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FINAL  REPORT 


on 

IMPLICATIONS  OF  ALLOY  STRUCTURE  WITH 
RESPECT  TO  MECHANICAL  PROPERTIES 
(Contract  N00014-76-C-0181) 

by 

L.  R.  Kahn,  A.  .1.  Markwortb,  P.  C.  Gchlen,  and  C.  T.  Hahn 
I . INTRODUCTION 

The  research  described  herein  was  carried  out  in  fulfillment 
of  Contract  N00014-70-C-0390  and  covers  rost-arch  carried  out  since  the 
last  major  rcport^'^  was  issued.  Publications  written  during;  this 
period  are  included  as  appendices. 

■'The  primary  objective  of  this  project  lias  been  to  relate  the 
toughness  of  metals,  including  effects  of  liydrogen,  to  physical  proper- 
ties at  the  atomic  level.  Our  approacli  is  based  on  computer -simu 1 at  ion 
techniques  that  describe  a crack-coutnin ing  crystal  comprised  of  the 

order  of  a thousand  atoms.  Witli  this  approach  (the  details  of  which 

(2) 

have  been  described  elsewhere  ),  each  atom  is  treated  as  a point  mass, 
and  appropriate  pairwise  central  potentials  arc  used  to  describe  the 
interaction  between  any  two  atoms.  The  discrete  atomic  array  is  inter- 
faced in  such  a manner  with  a surrounding  continuum  (see  Section  IV) 
tliat  its  behavior  is  characteristic  of  that  of  a corresponding  array 
contained  within  a macroscopic  crystal.  In  this  manner,  atomic  arrange- 
ments can  be  simulated  that  are  not  amenable  either  to  direct  experi- 
mental observation  or  to  treatment  via  closed-form  analytical  methods. 

A number  of  basic  problems  were  encountered  in  simulation 
studies  carried  out  under  a previous  contract^^^.  For  example, 
calculations  of  the  energy  barrier  to  crack  extension  incorrectly 
predicted  that  the  barrier  failed  to  disappear  at  large  levels  of 
applied  stress.  These  problems  were  ultimately  traced  to  two  fiindamcntal 
shortcomings  of  the  computational  procedure: 
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Enor^v  Bxchani;e.  The  proccdui'i'  failed  to  accoiml  for  the 
energy  exclianged  between  the  continmim  and  tlie  discrete 
atomic  array.  This  basically  was  the  reason  for  tlie 
anomalous  behavior  of  the  energy  barrier. 

Pi  sreuistry . We  liad  not  avoided  small  disregistr ies,  of 
the  order  of  the  atomic  spacing,  between  the  crack-tip 
position  implied  by  the  continuum  and  that  assigned  by 
the  discrete  atomic  array.  The  associated  errors  in 
energy  calculations  could  be  comparable  in  magnitude  to 
the  energy -exchange  discrepancy,  but  of  opposite  sign. 

In  some  cases  the  two  errors  canceled.  Wlien  they  failed 
to  cancel,  they  were  attributed  to  the  wrong  source. 


We  were  thus  led  to  believe  that  a small  rigid -boundary  model 
could  not  be  used  to  obtain  physically  meaningful  results;  effoi  s were 
consequently  applied  to  the  more  complex  flexible-boundary  model,  which 
contained  Che  same  flaw  but  for  which  its  diagnosis  was  more  difficult. 
During  the  present  reporting  per iod, methods  for  accounting  for  energy 
exchange  and  avoiding  disregislry  were  developed. 

Some  other  basic  questions  regarding  our  model  still  exist, 
e.g.,  concerning  the  manner  in  which  atoms  in  the  boundary  contribute 
to  the  overall  energy  of  the  discrete  region.  However,  wo  are  now 
confident  that  the  model  as  presently  constructed,  and  upon  which  the 
results  presented  lie  low  are  based,  constitutes  a useful  tool  for 
approximating  the  behavior  i^f  a "real"  material. 

"IXiring  this  reporting  period,  rcsearcli  was  carried  out  in 
the  following  areas:  (1)  construction  of  a piiysically  realistic  model 

for  the  atomic  arrangement  around  a crack  tip,  as  well  as  the  develop- 
ment of  "flexible  boundaries";  (2)  selection  of  an  Fe-li  interatomic 
potential  together  with  its  application  to  the  simulation  of  hydrogen- 
assisted  erttek  priipagat  ion ; and  (3)  analysi.s  of  the  effects  of  lattice 
periodicity  on  fr.acturc,  including  effects  due  to  the  presence  of  hydrogen. 
The  rescilts  are  presented  below. 


II . HYDROGEN -ASSISTED  CRACK  PROFACATION 


A.  Background 

The  presence  of  certain  foreign  species  (e.g.,  sodium,  helium, 
hydrogen,  etc.)  within  materials  can  result  in  effects  such  as  embrittle- 
ment, s tress "corros ion  cracking,  and  corrosion  fatigue,  which  are 
immensely  important.  One  aspect  of  this  very  broad  subject  that  was 
considered  in  this  program  was  hydrogen -induced  embrittlement.  Actually, 
as  has  been  pointed  out  by  Louthan  et  al.^'^\  hydrogen  embrittlement 
can  result  from  a number  of  processes,  including  the  following; 

1.  Hydride  formation 

2.  Reaction  between  liydrogen  and  some  impurity  or 
alloy  addition  contained  within  the  metal 

1.  Adsorption  of  hydrogen  onto  or  absorption  into  the 
metal 

4.  Hydrogen  "blistering"  or  cracking  (resulting  from 
hydrogen-gas  build-up  at  microcracks  and  voids). 

Mechanisms  proposed  to  explain  hydrogen  embrittlement  via  adsorption 
and/or  absorption  include  the  fo I lowing^ ; 

1.  Lowering  of  surface  energy  by  hydrogen  adsorption 

2.  Decrease  of  binding  energy  resulting  from  inter- 
action of  hydrogen  with  d-shell  electrons 

1.  Buildup  of  internal  gas  pressure  due  to  hydrogen 
precipitation  at  internal  voids. 

Absorbed  hydrogen  can  also  interact  with  dislocations  (e.g.,  acting  to 
pin  dislocations  or  loading  to  the  formation  of  gas-filled  voids  or 
microcracks  at  dislocation  pileups)  and  can  contribute  to  embrittlement 
in  this  manner. 

Our  studies  of  hydrogen  embrittlement  were  concerned  with 
effects  of  adsorbod/absorbed  hydrogen  on  crack  propagation  in  0(-lron. 

No  account  was  taken,  however,  of  effects  due  to  the  presence  of 


liyJrouen  on  the  Fc-Fc  interaction,  so  that  possible  contributions  due  to 

mechanisni  2 were  not  included.  Quantum-mechanical  approaches  are 
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currently  beiny,  used  at  Batte  1 le-Columbus  to  gain  information  of 
this  nature.  The  present  research  was  based  on  current  knowledge 
regarding  pairwise  Fe-Fe  and  Fe-H  interactions.  The  H-H  interaction 
was  not  included. 


B . Interatomic  Potentials  for  t'le  Iron-Hvdrogon  System 


Simulation  studios  of  hydrogen-assisted  crack  propagation  in 
<7-iron  were  carried  out  as  part  of  this  research  using  a simple  model, 
as  described  below.  As  far  as  the  Fe-H  interaction  was  concerned, 
initial  studics^^^  were  based  on  data  calculated  by  Walker  et  al.^^^ 
for  the  Fell  molecule  using  many-body  perturbation  theory.  I.ater 
studies  (summarized  below  and  described  in  chiail  in  Appendix  A,  and 

(B) 

Reference  (7))  were  based  on  the  Fe-ll  poienii.il  calculated  by  Olander'' 
using  exper  imen  t.i  1 1 V determined  data  for  hydri'gen  in  iron,  namely,  the 
heat  of  solution,  the  heat  of  .idsorpt  iv'n , .ind  the  tictivation  energy  for 
bulk  diffusion.  These  potentials  arc  presentei'  in  Figure  1 together 

(9) 

with  the  Johnson  Fe-Fe  potcnti.il  that  w.is  used  in  the  same  analyses. 
It  can  be  seen  tii.it  the  well  depth  for  the  Walker  potential  is  about  an 
order  of  magnitude  greater  th.tn  th.it  fi'r  the  Ol.inder  potential  and  also 
that  the  interatomic  distance  ;it  which  this  depth  occurs  is  much  shorter 
for  the  former  potential  than  for  the  l.itter.  In  addition,  the  Olander 
potential  extends  to  significantly  larger  interatomic  distances  than 
the  Walker  potential. 

It  must  be  emphasized  that  the  Fe-H  potential  is  of  great 
importance  in  this  type  of  study  .tnd  must  be  known  with  more  certainty 
than  is  presently  available  before  more  confidence  can  be  had  in  the 
results  of  the  simulations.  Nevertheless,  the  significant  differences 
between  these  two  Fe-ll  potentials  can  be  exploited  in  obtaining  infor- 
mation regarding  effects  of  variations  in  Fe-H  interaction  parameters 
(well  deptb,  etc.)  on  overall  crack-tip  behavior  under  the  influence  of 
applied  siress. 


Interatomic  Distance, 


FtCURE  1.  II.I.USTKATION  OK  INTKRATOMIC  POTKNTIAI.S  USR!)  IN  THIS  KESEAROM 
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C . S imu 1 a t ion  Studie s 

For  our  analyses  of  hydrogen-assisted  crack  propagation a 
two-dimensional  crack  model  for  a-iron  was  used,  the  crack  "tip" 
actually  being  an  infinitely  long  "line".  Our  results  were  found 
to  vary  only  slightly  with  the  type  of  boundary  condition  used,  i.e., 
flexible  or  rigid.  A single  atom  of  hydrogen  (actually,  a "line"  of 
hydrogen  atoms  parallel  to  tlie  line  defining  the  crack  tip)  was  placed 
within  the  lattice  near  tiie  crack  tip  in  such  a location  that  maximum 
effect  on  crack-tip  behavior  under  the  influence  of  applied  stress  would 
likely  be  observed.  In  tliis  light,  this  system  constitutes  a model  of  a 
form  of  embrittlement  resulting  from  "absorption"  of  hydrogen  into  tlie 
lattice.  It  is  clear,  liowever,  that  tliese  studies  represent  only  a first 
step  toward  the  solution  of  a very  complicated  problem. 

For  tile  case  ( Re  1 er enci.'  5)  in  which  the  Walker  Fe-H  potential 
was  used,  it  was  found  that  the  effect  of  tiie  hydrogen  was  to  shorten  the 
distances  between  the  four  iron  atoms  which  surrounded  it.  Consequently, 
the  Fe-Fe  bonds  beyond  were  elongated,  and  finally  broke,  one  at  a time, 
at  a load  about  10  percent  smaller  than  the  normal  failure  load  for 
u-iron.  This  behavior  can  be  understood  in  terms  of  the  forti  of  the  Fe-ll 
potential  function  (see  Figure  1)  which  exhibits  a deep  attractive  well 

O 

located  at  a relatively  small  interatomic  separation  (1.49A). 

For  the  case  in  which  the  Olander  potential  was  used  (Appendix  A, 
Reference  7),  tlie  results  were  quite  differenl  . The  iiydrogen  was  found  to 
cause  severe  distortion  of  tiie  lattice  in  its  vicinity,  the  crack  tip  being 
markedly  widened.  The  degree  of  crack-propagation  cniiancement  was  found, 
however,  to  be  relatively  insensitive  to  the  level  of  applied  stress,  in 
contrast  to  the  corresponding  results  obtained  using  the  Walker  potential. 
The  distorting  effect  of  the  Olander  potential  can  be  attributed  to  the 
fact  that  it  becomes  repulsive  at  relatively  large  interatomic  separations 
(see  Figure  1),  the  effect  of  the  hydrogen  thus  being  to  "push  away"  the 
surrounding  iron  atoms. 

It  is  interesting  to  observe  that  both  the  Fe-H  potentials  did 
result  in  effects  wliich  could  be  interpreted  as  forms  of  embrittlement 
despite  the  fact  that  the  potentials  themselves  had  quite  different 
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properties,  as  illustrated  in  Figure  1.  We  are  thus  led  to  the  con- 
clusion that  the  Fe-ll  interaction,  as  well  as  changes  in  the  Fe-Fe 
interaction  due  to  the  presence  of  hydrogen,  must  be  known  with  signi- 
ficantly greater  precision  before  more  definitive  conclusions  can  be 
drawn  regarding  mechanisms  of  hydrogen  embrittlement  at  the  atomistic 
leve 1 . 

The  hydrogen-embrittlement  problem  is  examined  from  another 
point  of  view  in  Section  III.  Analyzed  there  are  effects  of  hydrogen, 
situated  at  the  tip  of  a crack,  on  the  energetics  of  crack  propagation 
through  the  discrete  crystalline  lattice. 

III.  EFFECTS  OF  LATTICE  PERIODICITY  O.N  FRACTURE 
A . liackui'ound 

For  a crack-propagation  tlieory  based  purely  on  continuum 
mechanics,  energy  changes  associated  with  crack  extension  (or  healing) 
are  independent  of  the  structure  of  the  crystalline  lattice.  However, 
since  a "real"  solid  exliibits  atomic-level  periodicity,  one  would  conse- 
qtiently  expect  periodic  effects,  on  the  atomic  scale,  associated  witii 
crack  propagation.  Such  effects  are  well  known,  of  course,  for  other 
types  of  defects,  such  as  dislocations  and  vacancies,  and  indeed  have 
also  been  described  in  terms  of  effeits  (m  fracture  (as,  for  example, 
in  References  10  to  18). 

Interest  in  lattice-periodicity  effects  results  principally  from 
the  fact  that,  for  a range  of  stress  intensities,  an  activation-energy 
"barrier"  must  be  overcome  in  order  to  extend  or  close  a crack  by  a single 
lattice-repeat  distance,  as  illustrated  in  Figure  2a.  The  magnitudes  of 
sucli  barriers  would  be  dependent  upon  the  stress-intensity  level 
as  shown  schematically  in  Figure  2b.  IVo  successive  minimum-energy 
configurations  are  illustrated  in  I igure  3 for  a crack  tip  in  a model  bcc 
crystal.  The  crack-tip  separation  for  these  two  cases  is  one-half  the 
lattice  parameter.  It  is  apparent  that  crack  extension  can  thus  take  place, 
from  one  equilibrium  site  to  another,  as  a thermally  activated  process. 


1 


I/) 

c 

0) 


(J) 

c/> 

<D 


cn 


u 

< 

0£ 

CJ 

X 

H 


H H 
Z M 
W 3 
H 

O V) 
a-  cc: 
CjJ 
CjJ 

z Di 

h-t  0£ 
^ < 
.J  m 

H >• 
CO  O 
oi 
as  w 
'O  z 
w 

CJ  » 
::::  Z 
H O 


Uh  F— ' 
C CJ 

z 

C CJ 
»— I z 
H F- 
< 

H-f  Pl, 

> CJ 
Q 
Lj 

r:  H 
H t--i 
Z 

o 

O ^ 


O CJ 


c 

o 

H 

H 

■*— 

CJ 

/-s 

CO 

/-s 

H 

_c 

Q 

CO 

CL 

'w' 

:r> 

a 

CL 

hJ 

M 

I— 

CJ 

* 

z 

JsC 

H 

o 

o 

<; 

O 

H 

w 

CJ 

f— » 

o 

IS 

CO 

u 

C 

CO 

cu 

A5J9U3  |DI|U0iOcJ  |D|01 


CJ 

o 

I— J 

Uh 


9 


FIGURE  3.  ILLUSTRATION  OF  THE  ATOMIC  CONFIGURATIONS  AROTOD  A CRACK  TIP  (K  = 0.2eV/A 
FOR  TWO  DIFFERENT  CRACK-TIP  POSITIONS:  (a)  crack-tip  located  at  Position 
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analogous  to  other  types  of  defect  migration.  Also,  herein  lies  an 
additional  possible  mechanism  for  embrittlement  by  a foreign  atomic 
species,  namely,  through  the  action  of  such  a species  (a)  to  reduce  the 
activation-energy  barrier  for  crack  extension  and  thus  permit  the  crack  to 
advance  at  lower  stress-intensity  levels,  and  (b)  to  enhance  thermally- 
induced  extension.  This  could  take  place  by  effects  that  the  foreign 
species  has  on  metal-metal  bontls  and/or  by  alterations  produced  in  atomic 
arrangements  at  the  crack  tip  due  to  the  presence  of  the  foreign  species. 

B . Simulation  Studies 
1.  (h'Tipu  tat  Iona  1 Procedure 

Our  simulation  studies  of  effects  of  lattice  periodicity  on 
fracture  were  based  on  a two-dimensional  crack  model  similar  to  that 
used  to  simulate  hydrogen -embri tt lemen t effects  (see  Section  IIC) . 

However,  as  indicated  above,  additional  complications  were  involved  in 
the  latt ice -periodic i ty  studies  due  to  disregistry  and  energy-exchange 
problems  as  well  as  to  problems  involved  with  determination  of  saddle-point 
configurations. 

A relatively  simple  model  was  found  to  be  adequate  for  these 

computations;  it  consisted  of  165  movable  atoms  surrounded  by  a rigid 

boundary.  Within  the  boundary,  the  atomic  positions  were  taken,  for 

any  given  location  of  the  crack  tip,  to  be  in  accordance  with  the  dis- 

(19) 

placement  field  predicted  by  linear-elastic  theory  . The  potential 

energy  of  the  atomic  array  was  determined  as  a function  of  crack-tip 

position  by  reducing  the  discrete  configuration  to  its  minimum  potential 

energy  for  various  locations  of  the  crack  tip  along  a <10O>  direction. 

The  relaxation  procedure  was  carried  out  using  the  "kinet ic -energy  quench" 

(2) 

method  described  elsewiiere  ■.  (Wo  note,  however,  tiiat  other  appropriate 
minimization  techniques  could,  iu  principle,  be  used  to  find  minimum-energy 
con  f igura t ions . ) 

The  energy  of  the  entire  crystalline  system  was  estimated, 
for  any  given  crack-tip  position,  by  adding  to  the  potential  energy  of 

the  relaxed  discrete  region  the  appropriate  contribution  of  the  linear- 

2 

elastic  region,  that  being  vf(Kb/K^_)  where  K is  the  stress  intensity, 
y and  h are,  respectively,  tlie  surface  energy  and  lattice  parameter  of 


1 

J 


r 


the  iron  orysLil,  K is  tlio  stross  intonsiiy  aL  Lhc  Oriffith  lovi‘l 

„ 0 o 

(n  111110  r i ca  1 valiios  lor  y , h , aiul  K,  lioino,  (l.OiSloV/A  , 2.8(i  A,  and  0.2A2 

o5/2 

oV /A  , ros|H'o  I i VO  I o')  , and  i is  a iinitloss  paramoLor  oqiial  to  Uio  numbor 
of  half  1 at  t i 00 -paramo  t o rs  tlio  crack  Lip  lias  been  moved  rolative  Lo  some 
arbiLrarily  clioson  roforonco  position  (f  bcins  positive  for  crack  healing, 
negative  for  crack  extension). 

The  atomic  configurations  which  were  sought  were  those  which 
exhibited  a minimum  in  potential  energy  as  a function  of  crack-tip 
position,  since  it  was  these  configurations  for  which  no  disregistry 
existed  between  tlie  crack-tip  position  implied  by  the  continuum  and  that 
assigned  by  the  discrete  atomic  array.  The  nature  of  tlicsc  minima  was 
found  to  he  dependent  upon  stress  intensity,  that  is,  an  energy  minimum 
might  exist  either  for  the  discrete  region  or  for  the  total  system 
(iliscrete  region  plus  continuum)  depending  upon  the  magnitude  of  the 
K-lcvel.  This  effect  was  an  artifact  of  the  model. 

As  expectc'il,  ail  jaceiU  min  ini. i obtained  through  this  procedure 
were  found  to  be  characterized  by  crack-tip  positions  separated  bv 
about  b/2,  which  is  the  lattice  repeat  dist.ince  along  a • 1 00 d i rec  t i on 
for  thi.s  crack  conf  i giirat  ion  . In  otlicr  words,  periodic  variations  of 
energy  associated  with  the  cr.ick  exliibiied  the  periodicity  ol  the 
lattice  itself.  Dev'iations  from  this  hiduivior  could  exist,  liowever,  if 
the  lattice  were  to  be  disturlied  in  some  m, inner,  such  as  by  the  presence 
of  a foreign  species  (as  was  indeed  found  i or  the  situ.ition  described 
below  in  which  liydrogen  was  placed  near  the  crack  tip).  It  was  the 
atomic  con  f i g.u  ra  t ion  s .it  uljacent  minima  that  were  ol  concern  in  subse- 
quent calculations,  as  now  described. 

Poten t i 1 I -energy  values  lor  the  atomic  array,  at  given 
positions  along  .i  vector  in  con  i i gnr.i  t i on -.^pace"  |oining  two 
adjacent  minima,  were  ost  imatoi.1  by  using,  a procedure  in  which 
each  atom  in  I he  array  and  in  the  boundary  was  "marched"  (i.e.,  displaced) 
bv  identic.il  fractional  amounts  along  t lie  vector  (in  three-dimensional 
space)  Joining  its  positions  at  the  two  minima.  Tiu'  potential  energy  of 
tile  ariMv  was  calculated  at  tliese  various  intermediate  positions  along  the 
"march".  The  ma::  i mum  energy  thus  calciiiated,  at  one  of  tlie  intermediate 
positions,  represented  an  estimate  (actually,  an  overestimate)  of  the 
"sadd  le -po i n t " energy  of  the  discrete  array. 


For  a );i  veil  N-atom  system,  its  vector  in  configuration  space  would  con- 
sist of  the  set  of  3N  position  coordinates  cliaracter fz ing  all  the  atoms. 
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Two  points  shouKl  be  noted  relative  to  the  "march"  procedure; 

(11  Only  the  twi'  oiui-point  cim  fiyurat ions  of  atoms  were  "relaxed"  (i.e., 
minimum-ener,y\’'»  cu'ii  f i yura  t ion  s . (2)  One  can  siiow  tliat,  in  the  linear- 

elastic  approx ima t ion , the  fraction  which  describes  the  degree  to  whicli 
the  atoms  are  displaced  between  the  end  points  equals  the  fraction 
describing  the  degree  to  which  the  crack-tip  (as  defined  by  the  linear- 
elastic  field)  is  displaced  between  its  two  end-point  positions. 

Our  saddle-point  calculation  was  carried  out  using  an  approach 
based  on  the  technique  of  Sinclair  and  Fletcher  the  latter  having 

been  adapted  from  the  conjugate-gradients  metliod  for  function  minimization 
developed  earlier  by  Fletcher  and  Reeves^"  The  Sine lair-F letcher 
method  is  designed  to  locate  equilibrium  con  1 i gu ra t ion s ciiarac t e r i zed  by 
one  unstable  mode  of  vibration.  Two  attractive  features  of  the  method 
are:  (1)  Its  success  iloes  not  depend  on  constraints  placed  upon  any  of 

the  coordinates  or  upon  symmetry  within  the  structure.  (2)  First  (not 
secotnl)  derivatives  of  the  objective  function  are  required.  Fi'r  our 
applications,  the  initial  estimate  used  for  the  saddle-point  1. 'cation  was 
the  atomic  configuration  corresponding  to  thc'  maximum  ener,',y  I'f  the 
discrete  array  calculat:ed  using  the  "march"  procedure. 

2 . Results  - F.ffects  of  Vhariation  of  Stress  Intensity 


Computations  were  carried  out,  using  the  pri'ca'tlures  lieserihed 
above,  for  three  different  levels  of  stress  intensiiv;  K = 0.2,  O.J'ti.', 

o 5 /2 

and  0.4  cV/A  , the  second  of  the  three  values  bo  in;-,  the  stress-intensity 

level  K^,  prediett'd  from  the  Griffith  theory.  tine  vyould  expoi  t , on  the 

basis  of  continuum  theory  alone,  that  crack  extension  would  occ\ir  fv'r 

K > K^,,  and  crack  healing  for  K < K^.  Our  goal,  in  carrying  out  these 

computations , was  to  determine  the  lieight  of  ar t i va i ion -energy  barriers 

for  both  crick  opening  and  closure  as  the  crack  lip  was  moved  along  a 

straight  lim'  from  one  given  mi  n imum-enorgy  positic'n  to  the  next. 

o 5/2 

Results  for  K - 0.2  eV/A  are  presented  in  Figures  4 through 
b,  the  ati'mic  configurations  correspond ing  to  thc  two  calculated  minimum- 
energy  configurations  (i.e.,  at  crack-tip  locations  0 and  b/2)  being  in 
fact  those  shown  in  Figure  1.  In  particular,  .shown  in  Figure  4,  is  thc 
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CON’TI\'UUM)  WITH  CRACK-TIP  POSITIOX  FOR  K = 0.2c\7a5/2. 
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variation  of  potential  energy*  of  the  atomistic  region  witti  crack-tip 
position,  in  Figure  3 is  that  of  tlie  continuum,  and  in  Figure  6 is  that 
for  the  entire  crystal  (i.e.,  tlie  sum  of  tlie  two).  For  purposes  of 
comparison,  all  three  figures  are  drawn  to  the  same  scale  of  both  energy 
and  crack-tip  displacement.  A number  of  points  should  be  noted:  First, 

it  can  be  seen  from  Figure  4 that  tlie  minimum-energy  points  are  now 
shifted  away  somewhat  from  the  end-point  (and  relaxed)  configurations. 

Second,  the  average  energy  increase  associated  with  the  crystallite  in 
extending  the  crack  between  adjacent  minima  is  about  equal  to  the 
concomitant  increase  of  surface  energy,  indicating  that  the  model  is  well- 
behaved.  Third,  we  observe,  from  Figure  6,  that  energy  barriers  to  both 
opening  and  closure  do  indeed  exist.  Of  course,  there  is  a net  increase 
of  energy  as  the  crack  tip  advances  from  one  minimum-energy  configuration 
to  the  next,  since,  tlie  K level  is  below  the  Griffitli  value. 

Thus,  despite  tlie  fact  tliat,  at  this  K level,  crack  closure 
would,  on  the  basis  of  continuum  theory,  be  energetically  favored  to 
occur  at  anv  crack-tip  position,  tlie  present  model  exhibits  an  energy 
barrier  to  closure  of  about  0.05  eV . Accordingly,  crack  healing  can  be 
regarded  as  a thermally  activated  process,  with  a ''characteristic  temper- 
ature" being  defined  as  F/k  i00°C  where  E is  the  barrier  height  and 
k the  Boltzmann  constant. 

Results  analogous  to  those  presented  in  Figure  6 are  shown  in 

® 5 / 2 

Figures  7 and  8 for  K = 0.242  and  0.4  eV/,4  , the  former  corresponding 

to  the  Griffith  value.  For  both  situations,  energy  barriers  to  both 
crack  opening  and  closure  exist.  At  the  Griffith  level,  little  net 
energy  change  exists  between  adjacent  minima,  as  expected.  Also,  for 

o 5/2 

K = 0,4  eV/A  , the  net  energy  at  successive  minima  decreases  as  the  crack  tip 
advances,  also  as  expected.  The  fact  that  the  energy  minima  calculated 
at  the  Griffith  level  arc  not  at  the  same  value  of  energy  can  be 
attributed  largely  to  approximations  inherent  in  the  model.  For  example, 
some  portion  of  the  elastic  energy  attributed  to  the  contituium  was  con- 
tained within  the  discrete  reg.ion  and  was  thus  also  included  in  the 
[lotential  energv  cakul.ited  for  that  reg.ii'n. 

Absolute  energies  are  imimportant  tor  presiuit  considerations,  and  only 
energy  changes  relative  to  crack-tip  position  have  significance  here. 


0 0 


FIGURE  7.  VARIATION'  OF^ TOTAL  CRYSTALLINE  ENERGY  WITH  CRACK-TIP  POSITION  FOR 
K = 0.2Y2eV/A5/2_ 
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The  activation-energy  barriers  calculated  from  results 
illustrated  in  Figures  b tlirough  8 are  listed  in  Table  I.  The  barrier 
E_  for  backward  movement  is  seen  to  increase  rapidly  with  increasing 

D 

K-level,  whereas  that  for  forward  motion  decreases  with  increasing  K-level, 

albeit  more  slowly,  and  exliiliits  sonic  fluctuating  behavior  (E_  for  the 

o5/2  ■ °5/2  ^ 

K = 0.2  eV/A  being  greater  than  tliat  for  0.242  eV/A  ) due  to  the 

approximate  nature  of  the  model.  Also,  it  should  be  emphasized  that  these 

calculations  were  based  entirely  on  results  obtained  from  the  "march" 

procedure;  no  saddle-point  computations  were  involved  here,  although 

results  obtained  using  the  Sine lair-Fletcher procedure  were  obtained 

and  are  described  below. 


TABLE  I.  ACTIVATION -ENERGY  BARRIERS  CALCULATED 
FOR  F0RW<\RD  AND  BACKWARD  CRi\CK 
MOVEMENT 


05/2 

K(eV/A  ) 

Ep(eV)^^^ 

E^(eV)^^^ 

0.2 

0.  18 

0.05 

0.242 

0.27 

0.  '37 

0.4 

0.10 

1.54 

(a)  The  energies 
Figure  2a. 

E„  and  E„  arc 
r D 

defined  in 

Saddle-Point  Determination.  The  saddle-point  configuration 

could  not  be  determined  for  one  of  the  tliree  cases  studied,  namely,  that 
® 5 / 2 

for  K = 0.2  eV/A  since  the  crystallite  did  not  exhibit  a maximum- 
energy  configuration  at  some  point  between  adjacent  minima.  (A  maximum 
was  found,  however,  when  tlie  continuum-energy  contribution  was  added, 
as  illustrated  in  Figure  4,  so  that  the  absence  of  a maximum  for  the 
discrete  region  alone  was  simply  an  artifact  of  the  present  model.) 
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Nevertheless,  to  illustrate  the  procedure,  the  saddlc'point  configurations 

for  the  other  two  K-levels  were  deterniinetl.  Our  method  was  based  on  the 

Sinclair-Fletclior^'®^  technique,  and  included  modification  of  the  crack-tip 

position  at  various  places  in  the  saddle-point  search.  It  was  thus  found 

that  the  saddle  point  energies  for  the  discrete  region  for  K = 0.242  and 
o5/2 

0.4  eV/.\  were,  respectively,  only  about  0.02  and  0.06  eV  below  the 
corresponding  energy  maxima  determined  from  the  "march"  procedure,  so 
that  these  maxima  did  represent  relatively  good  approximations  to  the 
actual  saddle-point  energies. 

3 . Results  - Effects  due  to  Hydrogen  at  the  Crack  Tip 

The  computation  of  energy  versus  crack-tip  position  was  repeated 
for  K = 0.242  eVVA^*^^  but  witli  hyiirogen  present  at  the  crack  tip  in  a 
manner  analogous  to  that  in  an  earlier  study  (Reference  7,  Appendix  A). 

The  purpose  of  this  work  was  to  determine  the  effect  of  the  presence  of 
hydrogen  on  the  magnitudes  of  the  ac t i vat  ion -energy  barriers  for  crack 
extension  and  healing,  and  thus  to  assess  wliether  or  not  tlie  calculated 
alterations  due  to  hydrogen  would  result  in  effects  that  could  be  inter- 
preted as  embrittlement  of  the  metal. 

For  these  computations,  the  Olander  Fe-H  potential  ^ ' was  used, 
in  a manner  analogous  to  that  employed  in  previous  work^^^  . The  atomic 
configurations  around  the  crack  Lip  corresponding  to  two  adjacent  minima 
in  the  potential  energy  of  the  entire  system  (discrete  region,  including 
the  Fe-H  and  Fc"Fe  energy,  plus  continuum)  are  illusti'ated  in  Figure  9. 
These  minima  were  separated  by  slightlv  more  tlian  one  lattice  parameter, 
in  contrast  with  all  the  analogous  computations  without  hydrogen  for  which 
the  minima  were  separated  liv  only  about  half  a lattice  parameter,  an 
effect  brought  about  by  the  presence  of  hydrogen.  Also,  the  total  energy 
of  the  hydrogen -con t a i n i ng  system  was  found  to  bo  several  eV  less  that 
characterising  tlie  same  svstem,  at  the  same  K-lcv'cl,  but  without  the  "line" 
of  hydrogen  atoms. 

Results  of  carrying  out  tlie  "march"  procedure  are  presented  in 
Figure  10,  and  the  associated  barrier  heights  arc  listed  in  Table  II 
together  with  the  corresponding  energies  for  the  iron  system  at  the  same 
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FIGl'RE  9.  ATOMIC  CONFIGURATIONS  AROL^ND  A HYDROGEN -CONTAINING  CRACK  (K  = 0. 242ev/A5/2) 
CORRESPONDING  TO  ADJACENT  MINIMA  IN  THE  TOTAL  CRYSTALLINE  ENERGY. 
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K-level  but  without  liydrogon . It  can  be  seen  that  the  effect  of  the 
hvdroueti  is  to  decrease  and  to  increase  K,, . Consequently,  the  crack 

r 

is  energet ica  1 Iv  easier  to  open  and  more  difficult  to  cli’se  as  a result 
of  the  hydrogen,  factors  wliich  are  consistent  with  the  view  of  hydrogen 
acting  as  an  embrittling  agent. 


TABLE  II.  ACTIVATION -ENERGY  BARRIERS  AT 
K = 0.242  eV/A5/2  WITH  .\ND 
WITHOUT  THE  PRESENCE  OF 
HYDROGE.N 


Ep(eV) 

Eg(eV) 

No  Hydrogen 

0.27 

0.37 

With  Hydrogen 

0.16 

0.75 

.shouKi  be  noted  that 

tlie  specific 

results  obtained  here  are 

peculiar  to  this  particular  configuration  of  hydrogen  in  iron.  In 
general,  calculated  barrier  heights  are  expected  to  depend  upon  such 
factors  as  the  manner  in  which  hydrogen  is  distributed  along  the  crack- 
tip  surface  and  within  the  volume  of  the  iron  around  the  crack  tip. 
Additional  research  along  tliose  lines  is  needed. 

IV.  FLEXIBLE -BOUNDARY  SCHEME 

A problem  of  crucial  importance  in  the  atomistic  simulation 
of  crystal  defects  is  the  manner  in  which  the  discrete  atomic  array  is 
interfaced  witii  the  surrounding  continuum.  The  tise  of  "flexible 
boundaries",  in  which  atoms  within  the  bc’undarv  are  not  hold  rigidly 
in  place,  but  are  allowed  to  relax  tc'  positions  of  lower  potential  energy, 
results  in  a number  of  benefits.  Principally,  tlie  constraints  imposed 
by  "rigid"  boundaries,  which  can  leail  to  non-physical  beliavior  (such  as 
extraneous  stress  fields  at  the  boundary),  are  removed.  In  addition, 
with  flexible  boundaries  one  can  use  a smaller  array  of  atoms  in  the 
discrete  region  and  thus  improve  the  efficiency  of  the  computation. 
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In  t!\is  light,  an  updated  form  of  Battelle's  flexible- 
boundary  scheme,  "Flex-II",  was  developed  as  part  of  this  program,  and 
is  described  in  detail  in  Appendix  B (also  Reference  22).  The  method 
basically  consists  of  a set  of  steps,  for  most  of  which  certain 
alternatives  exist  regarding  detailed  procedure.  These  steps  are 
briefly  summarized  below  in  terms  of  a dislocation  simulation,  although 
the  procedure  is  applicable  to  other  types  of  line  defects  as  well. 

First,  an  array  of  atomic  positions  is  generated  for  a region 
of  perfect  lattice,  one  crystallographic  repeat  distance  thick  in  the 
direction  along  which  the  dislocation  is  to  lie.  Second,  an  initial 
dislocated  configuration  is  generated  by  displacing  all  atoms  according 
to  a linear-elastic  displacement  field.  T'ne  crystal  is  then  divided 
into  tliree  regions:  region  I containing  those  atoms  which  behave  in  a 

substantially  nonlinear  manner  and  are  thus  relaxed  independently; 
region  II  which  surrounds  region  I and  contains  all  atoms  which  lie 
witiiin  the  range  of  interaction  of  at  least  one  rcgion-1  atom;  region  III 
which  surrounds  region  II  and  is  comprised  of  that  portion  of  the  crystal 
which  is  necessary  to  fully  define  tlie  force  on  eacli  region-II  atom.  The 
following  third  and  fourth  steps  are  then  repeated  alternately:  Third, 

the  atoms  in  region  I are  relaxed  toward  cquilibri'im  according  to  the 
atomic-force  law.  Precise  equilibrium  is  not  sought  here,  since  later 
stops  disturb  the  region  bv  changing  the  botindary  conditions.  Fourth, 
the  forces  still  remaining  on  the  region-II  atoms  are  used  to  generate 
a corrective  displacement  field  for  all  three  regions  vising  Green- 
function  techniques.  The  last  two  steps  arc  repeated  iteratively  until 
no  atoms  in  regions  I or  II  have  a resultant  force  greater  than  some 
selected  tolerance. 


The  metliod  thus  devised  was  found 


(22) 


to  be  superior  to  earlier 


flexible-boundary  schemes  with  respect  to  speed  of  computation.  It  was 

(22) 

also  found  that,  in  the  analysis  of  core  contributions  to  a dislocation 
strain  field,  nonlinear  effects  were  fully  allowed  for  without  it  being 
necessary  to  use  a nonlinear  clastic  solution  for  the  displacements.  In 
addition,  this  approach  resulted  in  calculation  of  the  overall  dilatation 
of  a finite  body  due  to  a dislocation,  vjith  earlier  conflicting  calcula- 
tions of  this  quantity  liaving  been  corrected  and  reconciled. 
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V.  CONCLUSIONS 

• Effects  of  hydrogen  situated  near  the  tip  of  a crack  in 
O'-iron  have  been  investigated  using  two  widely  differing  Fe~H  pairwise 
potentials.  In  both  cases,  the  iron  was  ''embrittled"  as  a result  of  the 
presence  of  hydrogen,  although  the  mechanism  for  embrittlement  was 
different  for  the  two  cases. 

• Computational  procedures  have  been  developed  with  which  the 
energy  variations  associated  with  crack  propagation  through  a discrete 
crystalline  lattice  can  be  calculated.  These  were  applied  to  studies  of 
"lattice  trapping"  of  cracks  in  a~iron  as  a function  of  stress-intensity 
level,  including  calculation  of  ac tivation -energy  barriers  for  crack 
extension  and  licaling. 

• Effects  of  inclusion  of  liydrogcn  near  the  crack  tip  on 
lat t ice -trapp ing  properties  were  also  considered  for  a given  stress- 
intensity  level.  It  was  found  that  liydrogcn  caused  the  activation- 
energy  for  crack  extension  to  be  reduced. and  that  for  crack  healing  to 
be  increased,  factors  consistent  with  the  concept  oi  hydrogen  acting 
to  embrittle  iron. 

• A technique  was  developed  for  applying  flexible  boundaries 
to  atomistic  models  of  line  defects,  so  that  compatibility  can  bo 
achieved  between  the  discrete,  atomistic  region  and  the  surrounding, 
outer  region  described  in  terms  of  continuum  elastic  solutions.  The 
technique  was  found  to  he  superior  to  earlier  flexible-boundary  schemes, 
with  respect  to  speed  of  computation,  and  yielded  satisfactory  results 
relative  to  dislocation  properties  forOf-iron. 
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Ah St  rac  t 

Results  arc  presented  of  atomistic  calculations  of  the  eflects,  on  crack 
propagation  in  hody-contered-cubic  iron,  of  hydroj;en  atoms  located  at  the  crack  tip.  The 
Fe-Fc  interaction  was  described  in  terms  of  the  Johnson  inLeratomic  potential,  wliercas 
an  interatomic  potential  devclop^'J  ky  Olander  was  used  to  describe  the  Fe-ll  interaction. 
Tlie  11-H  interaction  was  not  considered. 

Calculations  were  carried  c'ut  using  both  rip.id  and  flexible  boundaries  and 
similar  results  were  obtained  for  the  two  cases.  A1  tl.oupji  the  presence  of  hydrogen  at 
the  crack  tip  did  indeed  :-esult  hi  an  enhancement  of  crack  prop.aga  t ion  relative  to  the 
situation  in  which  liydrop.m  was  n.>t  present,  the  degree  of  enhancement  was  found  to  be 
relatively  insensitive  to  the  level  of  externally  applied  stress.  This  contrasted  sharply 
with  the  results  of  earlier  calculations,  based  on  a different  Fe-II  potential,  in  which 
effo(ts  of  the  presence  of  hydrogen  were  found  to  vary  considerably  with  the  applied- 
stress  level.  The  differences  in  these  sets  of  results  are  indicative  of  the  existing 
need  for  accurate  interatomic  potentials. 

1 n t rodtic  t i C'n 

The  technological  importance  of  envi  i onmen  t a 1 1 y as.sisled  crack  propagation  is 
indeed  immense,  yet  tlu^  stale  of  otir  umler  st  and  i r.g. , at  a quantitative  level,  of  the 
phenomena  coitt  t i but  i ng  to  t li  i s effect,  is  relatively  primitive.  In  p.irticular,  the 
development  of  models  to  describe  the  detailed  n.iVure  of  the  discrete,  atomic  structure 
at  the  tip  of  .1  cr.uk,  whicli  is  under  the  influence  of  oxteinally  applied  stress  and  is 
peih.ii'S  anJuted  by  the  presence  of  f'  reii>,n  .itoms,  has.  become  a practical  possibility 
on  1 V witli  th*  advent  of  liirg.  -.sc.ile  lomputc'rs.  Ire  1 ini  nary  atomistic  calculations  have 
been  reported  elscwiu-te  (\)  c one  ern  iirg,  c'liecLs  of  hydrogen  on  crack  prop.igation  ina-iron. 
It  v.T.  thus  lovnrd  that  the  c i f ec  t -.f  ,,  hydrogen  atom,  situated  at  the  crack  tip,  was 
essentially  to  shorten  the  dislairces.  separa  I i np,  tlu'  four  iron  aloms  which  surrounded  it; 
consequen  1 1 V , oliier  Fe-Fe  h,.nus  were  elong.atc'd  and  tinally  broke,  one  at  a time,  at  a 
load  that  was  .ihont  10  percentt  l('Wer  than  the  norm.il  load  for  iailurc  of  rc-itrni.  A m.i  jor 
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Ehortccming  of  this  earlier  work  resulted  from  tlie  fact  that  the  Fc-H  interatomic  potential 
that  was  used  was  based  on  calculations  performed  fiir  the  Fell  mol  ecu  I e . The  research 
described  herein  consisted  of  additional  atomistic  calculations  of  the  effects  of  hydrogen 
o\i  crack  propagation  ina*ivon;  these  calculations  were  based,  however,  upon  a difiercnl 
Fe-1!  potential  t lian  the  one  nsid  in  the  earlier  work.  The  present  potcnli.il  characterizes 
the  pairwise  interaction  between  a hydrogen  .itoni  and  an  atom  of  iron  which  is  a con- 
stituent of  the  alpha  phase  of  the  metal. 

The  Fe-ll  I n t erne  t i on  i’ot  en  t i a 1 

Propcrticis  of  Fe -H  systems  have  been  the  subject  of  an  exten.sivc  amount  of 
research  and  have  thus  been  widely  discussed  in  the  scientific  literature.  Of  particular 
importance  for  atomistic  simulation  studios  is  the  interaction  potential  between  a 
hydrogeit  atom,  situated  within  the  iroit  lattice,  .'ind  any  given  iron  atom. 

Interaction  potentials  for  the  Fell  mo  1 ecu  1 e have  been  ascertained.  Thus,  for 
o.'camplc,  gas-phase  spec  troscopic  data  h.ive  been  useil  (2)  for  certain  of  the  first-row 
I rans i t i on-mo ta 1 hydrides  other  than  Fell  to  determine,  by  interpolation,  appropriate 
Morse-function  parameters  for  Fell.  .M  so , Knlker  ej_  aj_.  (3)  used  many-body  perturbation 
theory  to  execute  approximnie  lltir  t rci -lock  ca  1 cu  1 .i  t i i)n  s describing  a nundjcr  of  elc:Ctronic 
slates  of  the  Fell  molecule;  and  [ireviou.sly  described  ;;  imu  1 ;i  I.  i ons  (1)  of  the  effects  of 
hydrc'gcn  on  crtick  propngai  i o.-v  in  ot-iioi^  were  b.nsed  on  the  dat  a of  k'alker  c_l_  £2- 
the  ground-state  conf i gurat ion . 

The  solubility  of  hc’dro/.en  in  luilk  f ace -ecn  I erecl -c:u  hi  c metals  (including,  the 
gamma  pha;:c  of  irc'n)  has  been  detaribed  by  Vykbodet  s i._t^  aj^.  (4)  based  cm  a Morse-fenct  icin 
pairwit.e  potonlia.l  for  the  iiyd  rog.en -me  i a 1 interaction.  .A  Morse  pairwise  potential  lias 
al.so  been  used  by  Olandor  (3)  to  describe  hyd  rog.en -met  a 1 interactions.  In  Olandei  's  iv'ork , 
expo  r i men  t a 1 1 '.  (htermined  values  for  beat  of  solution,  heal  of  adsorption,  and  activation 
enc  >'g,y  for  bulk  diffusion,  all  for  hydrog.en  i ti  the  pertinent  metals,  were  tiscd  to  deter- 
mine' values  for  the  Morse-f  unc  t i on  p.i  i ,'ime  te  rs  for  the  various  liydrogen -mota  1 interactions. 
It  was  .in  Olander  Fc-ll  potential  which  was  used  in  tiie  present  study,  i.c.. 
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where  V(r)  is  tlu'  Fe-Il  inter. ici  ion  eneig.y  at  separation  distance  r,  with  r /a  - 1.12  (,a 
being,  the  l.ittice  parameter,  .about  2.86  A for  body-ccntcred-c  ubic  iron),  ra  = 2,96,  and 
1)  4.6  Kcal/ra'le  (0.20  cv)  . This  potenti.il  wa.s  determined  for  a { lOO]  surface  althougji 

v.iliiet:  for  these  pa  r.ime  1 1- r s were  fcniiul  to  vary  only  to  a miiv'r  dcg.rcc  with  choice  for  the 
ad.'.oi  fit  i on  plane.  For  our  cm  I c u 1 .i  t i on  s , the  pot  en  t i .i  1 was  set  equal  to  zero  for  Fe-H 

O 

sefi.i  ra  t ion  di  stances  greater  I b.an  about  20  A. 
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It  shoiiKI  l>o  noted  i li  U fundnmeitt.il  slvKlies  of  the  inlcrnction  betwctcii  hyilrojtt  n 
nnci  mot;!!',  are  t on  t i nu  i njt  niul  tho;:o  will  yii'ld  added  ;;  i j’o  i f i l ant  insiplils  int('  this  very 
ii»lH'rtanI  proMeni.  'Unis,  for  exa'np  1 e , Yin;',  cyt  aj_.  (6)  liave  develoiud  a density- 
lunelional  theory  I'l  >lier.u  sor  pt  i on  uainr,  a ;,e  I f -eons  i s t en  t.  1 i near -re  spon  se  formalism  and 
havi'  applii'd  their  analysis  to  the  ease  of  hydrogen  chemisorbed  on  a tunpsten  surface. 

Till'  interaeiion  enetp.y  between  the  hyd rc'p.iai  1 '"’i'  ‘>'"1  ll’*e  metal  surface  was  determined  as 
a function  of  distance  of  the  ion  iron  the  ;;nrfacc.  Likewise,  interatomic  forces  are 
heiret  calculated  by  Kahn  (2)  uainp  a noeel  theoretical  approach  based  on  recent  advances 
in  the  field  of  c|uantnm  chemistr\-.  Kahn's  initial  calculations  have  been  directed 
toward  de  t ermin  i ne,  interatomic  foiies  a.ssociated  with  a cluster  of  10  iron  atoms.  Plans 
for  rofininp,  this  qnan t nm-mechan i ca 1 approach  (2)  include  detailed  comparisons  between 
selected  calculable  and  exiio  r i men  t a 1 1 y measurable  bulk  properties  such  as  bull;  modulus, 
elastic  constants,  structure  factors,  and  Compton  profiles.  Calculations  of  interatomic 
forces  asr.ociated  Vi’ith  forcipn  atoms  (such  as  hydrogen)  contained  within  the  metal  will 
also  be  included. 

Compn  t a t i on  a 1 Procedure 

The  computer  a 1 r.or  i I h-.n  utilized  in  this  study  was  formulated  in  order  to  obtain 
solutions  of  the  Ntwlonian  equal  j. ins  ftir  a .set  of  N attmir  (considered  to  be  point  masses) 
in  the  im-edinte  vicinity  of  a delect.  It  is  usually  con\’onient  to  start  the  computation 
by  tl  i sp  1 c.c  ine,  the  atoms  from  their  perti'cl  lattice  positions  accordin;’,  to  the  equations 
of  linear  elasticity  for  the  p.iven  defect  (in  this  particular  ce.so  the  linear  anisotropic 
formalism  of  Sih  and  Liehowitz  (7)).  Once  the  atoms  are  so  positioned,  the  net  force  on 
each  atom  can  be  evaluated  aiul  insi'ited  Inti'  th.e  ecinalions  of  motion.  These  equations  aic 
then  solved  numerically  until  a stable  minimum  in  tin  potential  inergy  is  oht.ained,  or 
equivalently,  until  the  force  acting  on  every  one  of  the  M atoms  is  ;i  rlii  trar  i 1 y small.  A 
complete  description  of  the  algorithm  is  g,iven  in  Refcianci!  b . 

The  present  work  applies  to  a crack  length  of  2c  that  is  embedded  In  a large 
c'ontinnn;'i  to  which  a load  a is  applied  in  the  <-100>  e r vs  t .a  1 1 op,  raph  i c direction.  The 
crack  plane  is  a flOO]  ]->]anc  and  the  crack  'ini-  Is  iii'the  ■''.(,)0]>  direction. 

Since  we  restrict  ourselves  to  t hi'  case  of  plane  stiain  (i.e.,  all  displacements 
in  the  ::  direction  are  zero)  the  c r vs  t a 1 I og.ra  ph  i c periodicity  is  maintained  along  Llie 
crack  line.  Thus,  by  studying,  .slomic  displacements  in  only  one  e ry  s t a ] 1 og;r;iph  ic  repeal 
(listanee  eontaining  two  conseiutive  {OOlj  planes,  the  crystal  c.in  he  made  infinite  along 
the  crack  line  by  using,  periodic  homulary  conditions  v.liich  permit  s i mu  1 taneoiis  motion 
of  cqnivaleni  atoms. 

Ihe  ilisciete  reg.i  on  conl.ains  ahont,  7fK'  aliin;..  If  nvu  lo'scop  i c properties  of  a 
l)')dv  ,ir(  to  he  simula’ed,  1 h i •;  rrsion  must  he  tain'onnded  by  a miu  h larger  continnum  Viitli 
which  it  inter. iitfi.  The  InterK-ce  hetv/een  the  two  rogiona  is  refen'rod  to  as  a flexlbli; 
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bouiuliry.  In  I'sscnco,  the  f loxiblc'-bouiulnry  scluini',  FI^EX  II,  allows  the  continuum  to 
accomr.todat  (■  nonlinear  effects  arisinj’  in  llic  iirmieilialo  vicinity  of  the  crack  tip.  Tliis 
sclteme  i ilesciihoJ  hi' 1 ow . 

After  all  ac<'ms  have  been  d i sp  1 aceti  in  accord  witli  principles  of  linear  elasti- 
city, the  iH'dy  coniainiii!',  the  crack  line  is  divided  into  a discrete  atomic  region  (1),  a 
hi'undary  regii'n  (II),  and  a siirroundinp,  elastic  conliivium  (111).  The  boundaries  bet\7ccn 
the  three  rep.ic'ns  are  cylindiical,  tlie  axis  of  cacli  cylinder  being,  tlic  crack  line. 

The  essential  properties  of  t lie  three  rcg.ions  are  as  follows:  Region  I is 

comi'riscd  of  a discrete  array  of  atoms  tliat  contains  tile  defect's  centc'r.  Here  the  net 
intciatomic  forces  are  relaxed  by  int  eg,i  at  i ng,  the  erjiiations  c'f  motion  of  the  atoms.  Tlic 
atoms  in  Region  11  provide  the  link  between  Reg.ion  I and  tlie  continuum.  Here  net  forces 
.are  also,  calculated  a tomi  st  i c.i  1 ly  but  the  relaxation  of  these  forces  is  accomp  1 i siic  d by 
appropriate  linear  elastic  d i sp  1 accmi'ii  t fields.  The  continuum.  Region  III,  contains  only 
the  necessary  complement  of  atoms  for  calculating,  the  forces  in  Region  II. 

The  computation  scheme  proceeds  as  foJlovts:  atoms  only  in  Rog,ion  I are  relaxed 

for  a perii'd  of  time  approxiiiia  t e ly  equal  to  the  time  required  for  a sound  wave  to  travel 

from  the  crack  lino  to  the  interface  bt'tween  I and  II  and  hack.  Dviring  this  period  the 

boundaries  remain  fixed.  At  the  end  of  this  interval  th<’  net  forces  on  each  atom  in 

Region  If  arc  calculated  u.sing  the  interatomic  |)otential.  Because  periodic  boundaries 

arc  used  along  the  crack  line,  these  atomic  force::  can  be  con\’crt('d  directly  to  line 

forces  with  their  line  diri'ctions  parallel  to  the  crack  line.  With  these  forces  knov.’n , 

atoms  in  all  three  region:;  are  displaced  using;  the  li;iear  elastic  theory  described  in 

* 

liquations  5 and  6 of  Reference  9.  These  two  step:;  ;ire  repeated  alternately  until 
finally  tlic  forces  acting  on  tlic  atoms  in  hotli  Regions  1 and  11  liavo  hecotse  vanishing,ly 
small  at  which  point  final  ccpi  i 1 i hr  i urn  throiig.hout  the  entire  solid  extending,  to  infinity 
has  been  reached. 

Tlic  load  is  applied  to  the  discrete  region  hy  clioosing  tlie  appropriate  atomic 
displacements  for  the  starting',  linear  elastic  conf  i g,urat  ion . The  atomic  displacements 
are  proportional  to  the  stress  intensity  factor  fot  plane  strain,  Kj  , which  is  given  hy 

= cr,/?c  , 

where  c is  half  the  crack  length  and  CT  is  the  stress  applied  to  the  body.  If  a body  con- 
taining; a cT.^ck  is  (lerfectly  hritlle,  a statement  of  the  "tiriffith  principle"  is  that 


The  linear  el.istic  scheme  hroaks  down  in  the  i rimed  i ;i  ti’  vicinity  of  tlie  gmint  where  the 
force  is  applied.  The  d i ::p  I ■icemen  I ol  this  point  must  hi'  deterinineii  numerically. 
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crock  o::ten;;iou  occurs  wluii  Lho  sumnuH  [on  ot  Liic  c■|l.|^^',c  in  Llic-  clastic  strain  oiu'r;',y  in 
ihc  tiody  pU\s  the  work  done  on  the  body  duriii}’  crack  extension  is  equal  to  the  eneri^y 
required  to  form  the  now 
value  <'f  K 


(.irfaces  at  the  crack  tip.  'I'hc  corrcspoiidinp,  critical  Criifitti 

-J/2 

namelv  K,  , was  found  to  be  0.77  HNm 


Kesiil  t s 

The  compulations  described  herein  constitute  an  extension  of  earlier  work  (1) 
which  was  based  on  the  potcniial  developed  hy  Walker  e_L  aj . (t)  for  the  Fell  molecule . 

This  i-)Oiential  ( J)  was  characl  ei  i^ced  by  a minimum  eru  r^jy  value  of  19600  cm  ^ (i.e., 

2.43  ev  or  50.0  Keal/mole)  winch  was  found,  via  their  calculations  (3),  to  exist  at  an 
interatomic  separation  of  1.49  A.  Those  results  can  be  compared  to  tlie  Ol.mdcr  potential 
(5)  used  in  the-  present  work  for  wliicli  ihese  same,  two  paramet  rs  were  quite  different, 
i.e.,  0.20  ev  and  3.21  A,  respectively.  Shown  in  Fiipirc  1 are  results  of  a sample  calcu- 
lation from  the  earlier  work  showing,  the  lattice  distortion  caused  by  the  hydrogen  atom, 
i.e.,  at  irtiet  ing  the  4 iron  alom.s  around  it  and  Ihu.s  elongntitig  the  Fc“Fe  lionds  beyond. 
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Distortion  Around  fh'ack  I'ij)  Ha:, i d on  I’otenli  al  of 
W.ilker  ej^  aj  . ( '•)  for  the  Fe-H  I ii  t c r.tc  t ion  (from  Kefcrence  1) 
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As  can  lu'  soon  in  tills  fif’ui  c,  t >ui  hy >1  ro <itoi:i  i,,  closi'ly  siii  rounclL'd  hy  its  four 
nca ri-sL-ne i ^;hbor  iron  aLor.is.  The  lot  L- vaird  moLion  oL  atom  14  e:ausos  bond  LL  to  be 
sliolL-ht'd  and  eventually  to  break.  In  a similar  eonipui  ation  without  the  hydrofieu  present, 
bond  id.  is  not  brob.on  until  the  stress  on  the  crysuial  is  raised  to  Kj  = .'4.33  at  which 
point  unstable  crack  propasta t ion  occurs,  the  cleavage  pi  une  lu:ing  the  expe.rimcn tally 
obseiVed  llOOi  plane.  The  embrittling  effect  of  hydisjgen  l.s  not  evident  for  load  levels 
loss  than  = 1.83  The;;e  observations  substantially  ag.ree  with  experimental 

observation.  Subsequent  co.mjnit  at  ions , howevtir,  showed  that  the-  Walker  £L  .^A-  potential 
did  not  reproduce  many  of  t h>.'  known  bulk  properties  of  hydrogen  in  iron.  I'or  instance, 
it  jiredicts  a heat  of  solution  about  an  order  of  magnitude  larger  than  experimentally 
observed  and  a slightly  negative  molar  volume,  which  is  experimentally  known  to  be 
pcv'.itive.  In  an  effort  to  correct  these  Jiscrepanc  ie.: , Ulandcr's  pottaitial  was 
subs  1 i tti ted  for  Walker's.  Ihe.  former  had  been  matelii’d  (5)  to  the  heat  of  solution, 
resulting  in  its  shallower  wf-11  depth.  Although  ! ho  molar  volume  was  not  rc-calculated 
for  the  new  potential,  the  position  of  tlu.  miwir.:um  in  Ulandcr's  potential  curve  strongly 
indicates  that  the  molar  volume  will  le  positive.  Furthermoii , the  computer  algoritiim  was 
iiirproved  by  incorporating  a new  flexibie  boundary  schcjiie  vdiich  takes  into  account  ti'.e  fact 
that  the  crystal  is  a non-h.on.ogeneou,-.  (i.e,,  a cracked)  body.  In  the  previous  calculation, 
tb.e  linear  elastic  formalism  used  to  find  the  d i ;;pl  acemeu  t s associated  with  the  forcc.s  in 
Keg.iou  11  (see  above)  wa.s  for  a lu'mogenoous  body  and  the  scheme  used  to  make  tlie  crack 
St  rtuis- f ree  wa.;  only  approxiniate . In  the  present  scliemc,  this  ijoundary  condlticui  is 
satislicd  rigourously.  linally,  the  range  over  which  the  hydrogen  and  iron  atoms  interact 

O O 

was  extended  from  4 A in  the  old  model  to  about  20  A in  the  prestait  model. 

Kcsnlts  obtained  from  using,  tiie  new  computer  algorithm  and  tlie  Oiander  potential 
(.3)  are  presented  in  Figures  2 tlirougli  4,  all  of  which  repitsent  relaxed  atomic  configura- 
tion.^. In  particular,  illustrated  in  Figure  2 is  a crack  configuration,  with  no  hydrogen 
ptesent,  and  with  Kj  = 2.22  wiiereas  in  Figure  3 and  4 are  cracks  containing  hydrogen 

vitti  K = 2.22  Kj  and  1.48  k|,.<  rcsptc t ively . (\'"ipari;.on  of  Figures  2 and  3 shows  that 
1 •^0  '*0 

the  hydrogtui  causes  si’vere  distortion  of  tlie  lattice  in  its  vicinity,  tlie  position  of  the 
crack  tip  (as  defined  by  tlie  line  segments  joining  iron  atoms  at  the  surface  along  the 
crack)  being  nsarkcdl y widened  resulting  from  the  presence  of  tlie  hydrogen.  Interestingly, 
it  is  the  ri'pul'.  i-/e  portion  of  the  potential  which  dom.inatos  the  effect  in  this  work 
whereas  it  w.is  tlie  attractive  portion  which  dominated  the  effects  found  in  the  earlier 
work  ( 1 ) . 

Wtien  the  Oiander  potential  was  used,  ft  was  found  that  even  though  the  hydrogen 
• ili'r.i  was  ])lac's]  origin. illy  well  witiiiu  the  iron  lattice,  the  iron  lattice  appeared  to 
"reject"  the  h;.  Irogea  ator.j  .md  I'ositi'ui  it  out.side  its  surf. ice.  lu  the  original  Walker 
1 ee[ 'll  .11  ion,  the  opp  ' , i . ..luired,  n.imciy,  ilie  .st.irt  ing  po.-.itioii  of  the  hydrogo'ii  atom 
wa;.  jii.i  out  ide  til.  iro.i  l .ttice  an. I during  ivla.x.ition  the  .item  wa.s  moved  to  its  po.iition 
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FIG  A 

Effect  of  Uecrear.iug  the  Viiluo  to  1,A8 

shov.'n  in  Figure  1.  Originally  it  vaa  tliouglit  that  Lliis  rcjciction  of  the  hydi'ogcn  an 

artifact  of  tlu'  relaxation  t eclu:  Lquer , since  the  hottow  of  the  well  of  the  Fe-ll 

interaction  occurs  at  3.21  A (compared  to  an  t:quilibrium  di.'Uaiice  of  £ipiiro>;ii;-3tcly 

O 

1.A3  A in  ;in  uncrached  lattice),  extiL'isely  large  repulsive  forces  could  act  on  tlic  hydrogen 
.atom  during,  tlu  first  few  rclaxiition  step's,  wlien  the  .itoia  is  at  some  dist.anco  from  its 
eeju  i I i I'rium  jusition  witliJn  tlic  iron  lattice.  In  an  effort  to  correct  this  situation, 
•sever. il  altenipts  \,a're.  mad'-  to  dat'.pen  tlu-  atom's  motion  hy  increa.sing  its  mass  by  two 

ord.'rs  of  inagni  tiulo . This  dampen  in  j;  h.uJ  no  efftet  on  the  end  result,  indie, at  ing  that  t:hc 

ton  1 i j'.uiM  t ions  'dit'wn  in  Figuies  3 .and  A .ire  indeed  tin.'  true  ecjui  librium  coni  i gar.',  t ions 
eoi  1 e'.pond  i ng  to  the  Ulander  Fe-H  i nt  e laic  t iojt . Go;:’.i'.ari  ng  I'igures  3 and  A,  we  find  that 
tlie  'fieit  of  incre.isinr.  the  v.ilue  from  1.A8  to  2. .3,3  is  to  incro.ase  tli.'  cr.ich- 

t j j>  distortion  l’r(.iu;.-,li  t aliout  by  tiie  hydiogen  to  an  even  gre.ati'r  dep.ree,  .i:;  oiu'  w.'uld  liave 
('.-'.p  er  t ed . However,  no  I’e--le  bonds  weie  l>rolen  by  i ne  re.-e.  i ng.  the  value  to  this  liigber 
li'Vi'l.  Aloiip,  Mu'so  1 iiX's,  v-'e  noto  that  iu  the  results  o1  the  previous  v.otk  (1),  olonp.ated 
bonds  did  iu  bed  break,  one  at  a I lire,  a;i  the  v.alue  was  iucte.ised. 


Conclusions 


As  we  have  already  described,  the  effects  of  hydrogen  on  the  nature  of  crack 
propagation  in  a-iron  were  quite  different  in  the  earlier  work  (1),  compared  to  those 
presented  here.  In  the  former  case,  Fe-Fe  bonds  broke,  one  at  a Lime,  as  the  level  was 
increased.  In  the  latter  situation,  bond  breakage  occurred  upon  introduction  of  hydrogen 
into  the  lattice,  and  the  effect  of  increasing  the  level  was  not  to  break  additional 
Fe-Fe  bonds,  but  only  to  increase  the  lattice-distorting  effect  at  the  crack  tip  caused 
by  the  presence  of  the  hydrogen.  Unfortunately,  because  of  the  very  small,  atomistic 
scale  at  which  these  computations  were  carried  out,  we  have  no  direct  experimental  data 
relative  to  hydrogen-enhanced  crack  propagation  in  a -iron  that  can  lielp  us  determine  which 
of  the  two  sets  of  computations  was  the  mote  physically  realistic.  Thus,  the  need  for 
reliable  Fe-H  interatomic  potentials  must  indeed  be  placed  at  a high  level  of  priority. 

It  appears  that  theoretical  work  of  the  type  being  carried  out  of  Kahn  (2)  and  Ying  ct  al. 
(6)  represent  promising  means  for  reaching  this  goal. 

Of  course,  the  computations  described  here  and  in  the  earlier  work  (1)  constitute 
only  a beginning  of  atomistic-level  studies  of  environmentally-enhanced  crack  propagation. 
Other  effects  must  also  be  considered,  such  as  those  resulting  from  distributing  the 
hydrogen  all  along  the  crack  surface,  not  just  near  the  most  advanced  position  of  the 
crack  tip,  and  those  resulting  from  the  presence  of  otlicr  elements,  such  as  sodium,  that 
may  also  enhance  crack  propagation. 
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M cx  i i)  1 c_Bmuida  ry_  Coiulifious  and  Nonlinear  Gcoincti~ic  I'ffccts 
' ~in  Atomic  i s 1 ocal  i on  Mod  el  I i n;; 

by 

P.C.  Gclilen,  R.G.  Iloagland,  J.P.  llirth  and  J.E.  Sinclair 

A 1)  s t r a c t. 

A technique,  Plcx-Il,  is  described  for  aj)]ilying  flexible 
boundaries  to  an  atomic  region  in  computer  simulation  of  disloca- 
tions or  other  line  defects.  The  method  result.s  in  continuity 
of  c(\ui  1 i hr  ium , ui\der  the  clioscn  interatomic  potential,  across  tlic 
interface  between  tiic  atomic  region  and  the  outer  region 
described  in  terms  of  anisotropic  elastic  continuum  solutions. 

T!ie  tcchniciue  lias  high  numerical  efficiency.  It  is  shown  tliat 
when  the  crystal  is  initially  dislocated  according  to  the  Vol terra 
solution  for  di  splaccmejits , the  finite  strains  give  rise  to 
geometrical  I'.onlincar  effects,  usually  disregarded  in  linear 
elasticity,  wliicli  contribute  to  a volume  change  of  the  crystal. 
Allowance  for  this  effect,  and  for  elastic  nonlinearity  in  the 
crystal  beyond  the  boundary  region,  allow  tlie  overall  dilatation  of 
a finite  body  due  to  the  dislocation  to  be  rigorously  comp\ited. 

I-or  illustration  of  the  geometric  nonlinear  effect,  and  for 
comparison  with  earlier  modelling,  methods,  example  comjiut  at  i ons 
arc  given  for  the  [100]  edge  dislocation  in  a-iron. 
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I.  I nt  voduc  t ioi\ 


Tn  atomistic  computer  models  of  crystal  defects  sucJi  as 

dislocations  or  crack  tips,  incompatibilities  at  tlie  boundary 

of  the  atomic  region  will  arise  unless  a "flexible  boundary" 

f 1 - S') 

scheme  is  employed.  Several  sucli  metliods  have  recently 

been  developed,  and  liavc  improved  the  efficiency  of  the  calcula- 
tions by  permitting  the  use  of  smaller  atomic  arrays,  have 
yielded  descriptions  of  tlic  dilatational  and  higher-order  elastic 
fields  arising  from  nonlinearity  in  the  core  region,  and  have 
made  possible  calculations  of  defect  mobility  by  removing 
the  constraints  imposed  by  rigid  boundaries.  Various  modified 

versions  of  tlie  flexible  boundary  metliods  liave  been  briefly 
r 4 5 "I 

described'  ’ '.  Furthermore,  our  current  work  has  revealed  some 
subtleties  in  tlie  int  erjiretat  ion  of  the  core  contribution  to 
the  field  of  a dislocation,  connected  witli  tlie  finite  nature  of 
the  strains  in  tiie  initial  Vo  1 1 e r r a - so  lu  t ion  configuration. 

These  geometric  nonlinear  effects  lead  to  revision  of  the  earlier 
descriptions  of  dislocation  core  field  strengths,  although  tlie 
final  configurations  previously  determined  are  not  affected. 

In  view  of  these  changes  in  botli  method  and  interpretation, 
we  here  present  the  details  for  the  evolved  form  "Flex- 1 1",  which 
we  have  found  to  be  most  efficient  for  computer  defect  calculations, 
and  spell  out  the  steps  required  to  form  an  accurate  representation 
of  tlie  elastic  field  strengths,  including  volume  change, 
associated  witli  a dislocation  core.  A brief  comparative  account 
is  given  of  earlier  flexible  boundary  methods.  Some  selected 


results  a're  presented  to  jn'miit  com]Kirison  among  the  models. 


2. 


1 1 . Compu  t a t i oiia  1 F r p c c cl u r c 
A.  Flex-n 

The  method  consists  of  a number  of  steps,  for  most  of 
which  certain  alternatives  exist  as  to  detailed  procedure.  In 
some  cases,  tlie  clioice  is  a matter  of  taste  or  convenience. 

Somcwliat  surprisingly  liowevcr,  in  the  stop  of  sotting  up  initial 
positions,  wliat  appear  at  first  sight  to  be  arbitrary 
alternatives  lead  to  important  differences  in  effect.  For 
clarity,  we  first  give  the  basic  sequence  of  steps  and  then 
discuss  cac!i  step,  w i t li  any  alternative  options,  in  turn.  Although 
the  motliod  is  applicalile  to  dislocations,  crack  tips,  and 
otl\er  straig.ht  line  defects,  for  brevity  tlie  dcscri]>tion  is  given 
primarily  for  the  dislocation  case. 

The  first  step  is  to  generate  an  arra>'  of  atom  positions 
for  a block  of  ]U’rfect  lattice,  one  c ry  s t a 1 1 ogr  ajih  i c rcjieat 
distance  tliick  in  the  direction  (along  wliich  the  dislocation 
is  to  lie),  and  of  sufficient  size  in  tlie  x and  \n  directions 
to  accommodate  the  reg,ions  chosen  later  in  the  dislocated 
crystal.  This  stra  iglit  forward  step  will  not  be  dwelt  ujion  excejn 
to  remark  tliat  the  clioice  of  a higli  index  direction  for  x^  can 
result  in  the  formation  of  a series  of  jogs  or  kinks^^\  if  study 
of  these  is  di'S  i red  . 

The  next  stej)  is  to  form  an  initial  dislocated  configuration 
by  displacing  all  atoms  according  to  a 1 i near- c 1 a s t i c d isiilaccment 
field.  This  will  be  discussed  in  more  detail  below.  The 
resulting  array  is  tlien  divided  into  three  luirts,  as  in  Figure  1. 
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Region  I contains  those  atoms  (in  reality,  atom  rows)  which  arc 

' to  he  relaxed  i nclepeiulent  1 y because  tlieir  interactions,  according 

to  the  chosen  interatomic  iiotential,  are  substantially  nonlinear. 

Atoms  in  regions  IT  and  III  interact  via  tlic  same  potential, 

but  they  arc  always  displaced  collectively,  according  to  clastic 

Green  functions.  Region  II  contains  all  atoms  on  which  a force 

may  be  exerted  by  at  least  one  region- I atom,  while  region  III  is 

that  portion  of  the  remainder  of  tlie  crystal  wliich  is  necessary 

to  fully  define  tlie  force  on  eacli  region-II  atom.  Tl\c  thickness 

of  botli  regions  II  and  III  should  thus  equal  the  maximum  range 

of  the  interatomic  force  law.  Many  cmjnrical  potentials  for 

metals  liavc  been  f ormu  1 a t ed  with  a s])ccific  cut-off  at  less 

than  tlie  fourtli  nearest  neighbour  distance.  An  artificial  cut-off 

must  be  used  for  longer-ranged,  c.g.  ionic,  potentials'’^ 

The  shape  of  region  I may  conveniently  be  chosen  as  eitlier  a 

circular  cylinder  centerd  on  the  dislocation,  or  as  a roug!il>’ 

square  block  aligned  with  the  distorted  crystal  jilanes. 

The  following  final  two  stejis  are  repeated  alternately.  The 
first  is  to  relax  the  atoms  of  reg.ion  I tonuirds  their  equilibriur: 
state  according  to  the  atomic  force  law.  Since  later  stcjis 
will  disturb  the  region  by  changing  the  boundary  conditions,  this 
relaxation  stage  is  not,  for  efficiency's  sake,  greatly  prolonged 
in  search  of  precise  equilibrium.  If  a modest  limit  (e.g.  50) 
is  set  for  this  ste]i  on  tlie  numlier  of  force  calculations  allowed, 
then  in  early  iterations  a useful  economy  is  obtainctl,  while 
in  the  final  iteration  a sufficiently  precise  equilibrium,  results, 

I s ince  t he  f i na 1 relaxation  is  only  a small  adjustment. 
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The  second  of  tlic  iterated  steps  is  to  use  the  forces  still  ] 

remaining  on  region  -II  atoms  to  generate  a corrective  displacement 
field  for  all  three  regions  (and,  implicitly  , the  entire 
crystal),  composed  f rom  a suiierpos  i t ion  of  Green  functions.  In 
this  two-dimensional  context,  a Green  function  for  a lattice 
or  continuum  is  a displacement  field  wliich  will  produce  ccjuilibrium 
in  res]-)onse  to  a line  of  force  of  unit  strcngtli  in  a cliosen 
direction,  applied  at  some  origin.  (Alternately,  the  same 
displacements  remove  a line  of  force  existing  through  internal 
interactions).  It  is  thus  necessary  in  tins  step  of  tlie  calculation 
to  form  a sum  of  Green  function  fields  w'itli  origins  at  each 
of  the  region-II  atoms,  and  source  strengths  given  by  the  forces 
on  those  atoms,  first  converted  to  units  of  force  per  unit  x^" 
length.  Because  region  I is  nonlinear  and  inhomogeneous,  this 
step,  which  uses  the  Green  function  for  an  infinite  homogeneous 
body,  cannot  prod\ice  exact  equilibrium.  Iteration  of  the  final 
two  steps  is  tluis  necessary.  However,  to  the  extent  that  regions  II 
and  III  behave  linearly,  no  forces  requiring  correction  by 
Green  functions  should  ever  develop  in  region  III  or  beyond,  for 
all  the  interacting  neighbours  of  any  atom  in  region  III  are 
displaced  only  according  to  equilibiium  1 incar- e 1 as t i c fields. 

Tlie  two  final  steps  should  be  iterated  until  no  atom  in 
regions  1 or  II  has  a resultant  force  greater  than  some  cliosen 
tolerance.  Convergence  in  3 to  6 iterations  is  usual  in  our 
ex])erience.  During  tlie  iterations,  it  is  useful  to  accumulate 
certain  information  on  the  source  strengtlis  of  the  Green  function 
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displacements  ajiplied,  for  use  in  analysis  of  tlic  dislocation 
core  field.  This  will  be  discussed  in  section  ITB,  after  tlie 
following  enlargements  on  the  alternative  options  available  at 
eacli  step. 

I n i t i a 1 jio  s i t i on  s from linear  elasticity 

Regarding  the  choice  between  isotropic-  and  an  i so  t roji  i c - 

elasticity  solutions  for  the  field  of  a dislocation,  we  believe 

that  the  latter  results  should  be  used  in  all  cases.  Although 

trials  liave  shown  little  effect  on  core  structure  from  using 

tlic  inappropriate,  but  analytically  simpler,  isotropic  expressions, 

the  anisotropic  results  are  now  well  documented,  and  new,  num.crically 

efficient  methods  of  cominiting  tliem  liave  become  available. 

Arbitrary  orientations  and  anisotropy  can  most  expediently  be 

dealt  with  using  the  so-called  integral  f orma 1 i sm ^ ^ . Alternatively, 

(12  1?)') 

for  li i gli- symmetry  cases,  the  "sextic"  formalism  ’ using 
complex  functions,  leads  to  sim]ilc  expressions  wliich  may  be  easier 
to  implement  if  great  generality  is  not  required.  All  of  tlicse 
solutions,  whether  isotropic  or  anisotrojiic , we  shall  loosely  term 
"Volterra  fields",  meaning  that  only  the  dominant,  logaritlimic 
singularity  is  included. 

A second  ([uestion  whicli  has  received  discussion  in  previous 
dislocation  modelling  work,  is  the  possible  unexpected  asymmetry 
which  can  occur  wlicn  the  Volterra  solution  for  an  edge  dislocation 
is  aiiiil  ied  in  tlie  most  obvious  way.  Tliis  can  arise  because 
the  linear  elasticity  tlieory,  being  based  on  an  infinitesimal 
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strain  measure,  does  not  distinguisli  the  geometry  of  the  initial 
and  distortet!  states.  Thus,  the  loss  of  symmetry  which  occurs 
for  edge  dislocations  wlien  the  Volterra  solution  is  taken  as 
a function  of  initial,  or  Lagrangian  coordinates,  is  not  "noticed" 
by  the  infinitesimal  theory.  It  happens  that  this  asymmetry 
can  be  avoided  if  the  branch  cut,  or  surface  of  discontinuity  in 
the  I'olterra  field,  is  moved  from  the  axis  to  the  positive 

axis,  corresponding  to  dislocation  formation  by  insertion  of  a 
half-plane  of  atoms,  rather  than  by  slip.  This  was  the  approach 
used  by  Gehlen  et  al^^^.  An  alternative  sclieme^^^  which  results 
in  the  expected  symmetry,  is  to  use  final-state,  or  Hulerian 
coordinates.  To  do  tliis,  it  is  simply  necessary  to  iterate 
the  calculation  of  tlie  di  S]-)!  acement  of  each  atom  until  its 
displacement  from  tlie  corresiionding  {lerfect  lattice  jiosition  is 
sc  1 f- cons  i s tent  1 y given  by  tlic  Volterra  solution  evaluated  at 
the  displaced  position. 

The  clioicc  from  among  the  above,  or  otlier  similar  ojicrational 
schemes  for  setting  uji  initial  po  s i t ions^  would  appear  at 
first  sight  to  be  arbitrary,  since  all  are  equally  valid  within 
the  terms  of  the  i n f i n i t es i ma 1 - s t ra i n , linear  elasticity  theory. 
However,  wiicn  calculating  the  volume  change  associated  with  the 
core  field  whicli  is  found  in  tlie  covsrse  of  flexible  boundary 
relaxation,  it  turns  out  that  the  small  differences  between  the 
various  above  interpretations  of  the  Volterra  solution  cannot  be 
neg,  looted  , since  tliey  include  a significant  relative  difference 
in  dilatation,  both  local  and  overall.  This  is  of  second  order 
in  the  strain  magnitude;  that  is,  of  geometrical  non-linear  origin. 
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Tho  well  knoKii  result,  that  the  Vol  terra  solution  causes  no  | 

i 

volume  change,  is  true  only  to  lowest  order  in  the  strain 
magnitude.  The  results  of  choosing  different  representations 
of  the  Volterra  solution  will  be  demonstrated  in  section  IIIB. 

Whatever  the  choice,  Tlcx-II  produces  an  accurate  equilibrium 

configuration.  However,  because  of  the  greater  case  in  calculating. 

the  overall  volume  change  (section  I IB),  we  prefer  to  use  ' 

the  fulcrian  coordinate  aiqiroacli  to  calculating  initial  positions. 

Atomic  Relaxation  i 

We  have  used  two  different  methods.  One,  the  so-called 
dynamic  metliod^^^^,  is  well  known  and  is  liascd  on  an  integration 
of  the  equations  of  motion,  with  periodic  or  continuous  artificial 
quenching  of  tlic  kinetic  energy.  The  other,  known  as  conjugate 
gradients  ^ , is  less  mechanical  in  ph  i losopliy  , and  consists 

of  a scries  of  sea  relies  of  configuration  space,  systematically 
reducing  components  of  the  energy  gradient  vector  to  zero.  One 
or  other  of  these  methods  may  be  more  efficient  for  a given  initial 
configuration  and  jiotcntial,  but  the  difference  in  ])crformancc 
is  never  very  great. 

The  Cl rcc n Tunc  t i on 

Because  of  tlie  gross  distortions  in  tlic  core,  the  body  is 
elastically  inliomogcncous  , and  the  true  Green  function  depends  on 
both  source  and  field  position.  N'eglect  of  this  fact  allows  us 
more  simply  to  use  a function  only  of  the  relative  source-field 


position.  But  because  of  the  dislocation,  tlie  relative  positioning. 

of  two  atoms  cannot  be  uniquely  defined  in  terms  of  lattice 

indices.  It  is  therefore  necessary,  and  in  any  case  much  simpler, 

to  use  a continuum  Green  function  for  moderate  and  large 

separations.  The  anisotropic  elastic  displacement  field  of  a line 

force  is  given  in  the  integral  formalism  by  Asaro  et  or 

in  the  sextic  formalism  by  Stroli^^^\  Tor  higli  symmetry  cases^^  ^ 

f 21 

or  tlie  isotrojiic  approximation^  • simjiler  explicit  expressions 
arc  available.  All  of  tliese  solutions  diverge  logar  itlimical  ly  at 
the  line  of  force,  so  that  tliey  certainly  cannot  give  ttic  dis- 
placement of  tlic  row  of  atoms  on  which  foices  act,  nor,  probably, 
that  of  otiicr  closely  neighbouring  rows.  However,  the  no  i ghbouriiood 
of  any  region  -II  atom  is  essentially  perfect  lattice,  so  that 
a lattice  Green  function  may  be  used  for  near  neiglibours.  This 
can  easily  be  calculated  by  performing  a relaxation  of  a small 
block  of  perfect  crystal,  with  the  central  row  of  atoms  bearing 
an  externally  ap]ilied  force,  and  tlie  boundary  atoms  being  located 
according  to  tlie  continuum  Green  function.  The  resulting  dis- 
placements per  unit  applied  force  arc  stored  in  an  array  for  use 
in  I-'lex-lI.  Tlie  size  of  the  atomic  region  whicli  has  to  be  relaxed 
in  this  sub.sidiary  experiment  in  order  to  give  good  equilibrium 
will  indicate  the  range  at  which  the  continuum  Green  function 
must  be  replaced  liy  the  lattice  Green  function.  In  model  a-iron, 
a range  of  about  7>h fl  was  found  sufficient.  In  early  work  on 
i'l  ex- I 1 the  need  for  a lattice  Green  function  was  neglected, 

zero  value  being  associated  with  the  row  at  tlie  origin.  Slow 
convergence  was  still  achieved  in  tlie  ricx-II  procedure,  but  this 
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could  not  bo  guaranteed  in  general  of  course. 

An  alternative  to  ]irc-calculat  ing  tlie  near-neighbour 
lattice  Green  function  values  is  as  follows.  Only  regions  I 
and  III  are  displaced  in  response  to  forces  in  region  II,  using 
the  continuum  Green  function.  The  required  displacements  for 
region  II  are  then  obtained  by  relaxing  tlic  atoms  in  that 
region  holding  regions  I and  III  fixed.  This  procedure  was 
devised  to  deal  with  applications  to  cracl.s,  where  the  Green 
function,  including  nea r - nc i ghbou r lattice  values,  is  a function 
of  source  position  relative  to  the  crack  faces. 

I? . A n a 1 )•  s i s of  d i s 1 o c a t i o n c ore  f i e 1 c’ and  voltime  cliange 

At  tlic  comiiletion  of  tlie  Plex-II  procedure,  tlie  d i sj)l  aceinen  t 
field  of  the  entire  crystal  excluding  region  I is  given  by  the 
Vo  1 terra  field  plus  an  additional  field,  i elated  to  tlie  core 
conditions,  comprising  a sum  of  Green  funttion  terms.  This  core 
field  can  be  completely  described  in  terms  of  t!ie  total 
distribution  of  line  forces  whicli  generate  it.  Tlie  projierties 
of  this  distribution  whicli  relate  to  the  clastic  far  field 
arc  the  moments  of  the  distribution  about  tlic  origin.  Of  greatest 
significance  are  the  dijiole  moments 

Mi  = r.  I- . (n)  X (n)  , i , J = 1 or  2 (1 ) 

where  F(n)  is  the  force  per  unit  x--length  at  x(n),  the  n-th  atom- 
row  iiosition.  If  (X2>::^)  is  a mirror  plane,  only  M^^^  and  M^-,  will 
be  non-cero.  The  dipole  moments  can  be  related  to  the  dilatation 
of  the  body  which  the  core-field  terms  contribute;  higher  order 


iiiomcnls  r.ivc  further  inlornation  such  as  llie  positioning  of 
the  a]>parent  core-fieUl  origin  relative  to  tlie  dislocation 
line  (see  sect.  II  D). 

Unfortunately,  as  exiilaincd  in  the  previous  section,  the 
definition  of  tlie  Volterra  field  is  somewhat  arbitrary,  so 
that  the  above  characterization  of  the  additional  core-related 
field  in  terms  of  linear  elastic  solutions  is  not  unique, 
but  may  only  be  given  with  respect  to  a particular  representation 
of  the  Volterra  field.  Nevertheless,  it  is  possible  to  determine 
uni([ucly  such  overall  features  as  tiie  formation  volume  of 
tlie  dislocation,  as  wo  now  describe.  To  do  so,  we  must  take 
into  account  not  only  the  arbitrary  choices  made  in  constructing 
the  Volterra  field,  sucli  as  origin  and  branch  cut  positions, 
but  also  both  geometric  and  elastic  non-linearities.  However, 
it  is  not  necessary  to  derive,  nor  to  use  in  the  rel  a.xat  ion , 
a complete  non  - 1 i near - e 1 as t ic  solution  for  the  field  of  a dis- 
loca t i on . 

We  may  list  the  far-field  terms  contributing  to  the  overall 
dilatation  as: 

1)  the  initial  Volterra  displacement  field; 

2)  the  Green- function  terms  applied  to  regions  II  and  III  during 
relaxat ion ; 

7>)  additional  d i sp  1 acemen  t s which  would  be  needed  to  achieve 
exact  ecpii 1 i hr ium  throughout  the  material  beyond  region  II; 
these  are  necessary  because  the  nonlinear  terms  in  tlie  stress- 
strain  law  lead  to  both  a body  force  distribution  throughout 
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this  outermost  region,  and  to  outer  surface  tractions,  when 
the  strains  are  given  by  linear  clastic  fields; 

4)  any  misfit  at  the  branch  cut  (slip  plane)  in  the  initial 
Vol terra  field. 

Contribution  (4)  arises  only  when  Lagrangian  coordinates  arc 
used.  A slight  overlap  of  material  results  from  failure  of 
the  I.agrangian  Voltcrra  solution  to  satisfy  exactly  the  reepn  red 
displacement  jump  conditions  on  the  slip  plane.  This  oddity 
is  scarcely  noticed  in  the  atomic  modelling,  but  must  be  alloKcd 
for  in  order  to  obtain  the  same  final  result  independent  of 
the  initial  configuration.  To  avoid  tliis  complication,  we 
restrict  tlie  follov.'ing  analysis  to  T.ulerian  coordinates,  and 
give  expressions  for  the  three  contributions  to  the  overall 
volume  change. 

Firstly,  in  order  to  compare  the  volume  of  a body  with 
tlie  Fulcrian  Volterra  displacements  to  that  of  the  same  quantity 
of  material  uitliout  dislocation,  care  maist  lie  taken  to 
define  the  position  of  the  slip  surface,  and  fur  t licrmore , the 
following  exact  expression  for  tiic  area  cliange  of  a contour 
must  be  employed: 

du 

6A  = / I u x d)c  I - i / |u  X — |dJ,  (2) 

a ~ “a'  dz~ 

in  wliic!)  tlie  first  term  is  the  customary  infinitesimal  strain 
result,  and  tlie  second  term  enters  because  the  material  points 
connected  by  dl  in  tlie  I'inal  state  have  a different  se]iaration 
in  the  initial  state.  In  cq . (2),  I is  a closed  contour  in  the 


A 
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final  state,  and  u is  a displacement  field  expressed  in  terms 
of  final  coordinates.  When  ajiplied  to  the  Volterra  field  for 
an  edp.e  dislocation  in  a circular  cylindrical  body  of  radius  R, 
equation  (2)  yields  for  tlie  volume  change  per  unit  lengtli 

[a^2,n(R/r^)  + o^]  (3) 

Here  b is  the  Burgers  vector,  taken  to  be  in  the  direction, 

and  n,  and  ot , depend  on  the  elastic  constants.  For  isotropic 
e 1 a s t i c i t >• , . 

l-2v  . _ 5-4v 

“l  ■ TtiCI-vI  ’ “2  ■ SttCT^J  ' 

where  v is  Poisson's  ratio.  For  anisotroj-))' , is  simply  the 
coefficient  of  bJ.nr  ajjpearing  in  the  x ^ - d i sp]  acernen  t component; 
no  expression  for  «2  given,  since  it  depends  upon  the 

form  of  the  relaxation  terms  arising  at  a free  cylindrical  surface 
]. resent  ly  unknown  for  anisotropy.  Tlie  choice  of  slip  surface 
position  enters  through  r^ , in  eq.  [3),  defined  to  be  tlie 
distance  from  the  origin  at  which  the  sli])  surface  in  tlie  final 
state  crosses  the  x^  coordinate  axis.  This  may  be  altered,  for 
instance,  by  shifting  the  atomic  structure  relative  to  the 
coordinate  system  in  such  a way  that  climb  aiipears  to  take  place 
without  in  fact  altering  the  amount  of  material  at  the  core. 

Secondly,  the  volume  change  during  relaxation  by  Flex-II, 


which  we  shall  express  as 
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is  {;ivcn  in  terms  of  the  dipole  moments  (eq.  1)  by^ 


6A„,.  = S . . . , M. , 
Cl-  iijkjk 


C6) 


where 


■^ijkH  ^ '\jke,  ■ ^ij33^3k£'^‘‘’3333’ 


(7) 


and  ‘'inisotropic  compliance  constants  approjiriate 

to  the  clioscn  dislocation  orientation.  It  is  important  to 
realise  tliat  tliis  formula  applies  to  a finite  liody.  Tlie  "infinite 
body"  dilatation,  i.e.  that  of  a larye  material  surface  far 
from  botli  tlie  defect  and  tlie  free  surface,  is  related  to  tlu'  M^j 
by  different  compliances. 

Finally,  tlie  volume  chanye  associated  witli  liyimt  bet  ical 

relaxation  of  tiu’  nonlinear-clastic  effects  licyond  regions  I 

and  II  can  be  found  to  sufficient  accuracy  by  c.ilculating  the 

line  a r clastic  relaxation  of  the  remaining  forces  associated 

with  nonlinear  effects.  In  terms  of  the  third  order  constants 

1191 

d.  defined  bv  Willis  ' , tlie  body  force  per  unit  volume, 

cxistini;  under  Failerian  displacements  u^^  satisfying  the  linear-clastic 
equilibrium  equation,  is 


fi  = cK: 


9 |-9U(i  9u,,, 


ijkX,pi,i  9x^  9x], 

In  terms  of  moments  of  tins  distribution  defined  by 


(8) 


'^~rh~i"s  T‘i  T iu.~  re's  ull^f  or  j)  1 alie~srraln  and  anisotroiiy  corrcsiiond  ing 
to  the  three  dimensional,  isotrojne  result  given  on  page  89 
o f r.she  1 by  (ref.  18). 
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tlio  (' i n i t 0 - body  volume  eliangc  associated  witli  1 j nca  r - el  ast  ic 
relaxation  of  tliese  forces  is  given  by  the  same  formula  as 
Cl].  ((>)  • i >1  addition,  the  non-linear  clastic  law  will  lead  to 
surface  tractions  wliose  relaxation  would  iiroduce  a furtlicr 
dilatation.  Because  we  do  not  liavcjO\'en  for  tlie  simple  case 
of  a circular  cylindrical  body,  exact  cx])ressions  for  the 
strains  at  and  near  tJie  free  outer  surface  (if  tlie  elasticity 
is  anisotropic),  tlic  last  mentioned  contr  i iuit  ion  to  t!ie 
dilatation,  and  also  j-art  of  the  contribution  from  the  body 
force  (8),  c.mnot  be  determined.  However,  tlie  principal  contri- 
bution from  tlie  nonlinear  effects  can  l)o  found  by  using  the 
infinite  body  expression  for  the  strains  in  (8),  wliicli 

tabes  tiie  form 


,^'il  = In  ('01 


(10) 


wliere  r,  0 are  polar  coordinates.  Substitution  in  (S)  leads  to 


-.8 


fj(t,4^)  = r 'g-Cfl) 


(li) 


and  finally,  t!n'>uj;h  (0),  to  the  following  form  for  tiie  total 
non-linear  elastic  contribution  to  tlie  dilatation: 


" ^‘'(^.I'nUl/rj,)  + a^)  , (12) 

where  r^  is  tlie  inner  radius  at  wliicli  the  integral  (9)  commences. 
Section  111-C  slioiv's  liow  the  effect  iv'c  value  of  r^ 


c a n b e 
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determined  in  practice,  hut  clearly,  it  corresponds  closely 
to  tlic  inner  radius  of  region  III.  The  value  of  is  well 
determined  tlu'ough.  t!ic  second  and  third  order  clastic  constants 
entering  in  (8)  and  (10).  The  unknown  surface-relaxat ion 
effects  arc  emhraced  hy  tlic  term  a^.  Consideration  of  the  order 
in  r of  the  surface  terms,  and  examination  of  t!ie  isotropic 
case,  show  that  is  independent  of  R.  Hence  this  term  could 
alternatively  he  emhraced  in  the  rt^  term  by  taking  an  "effective" 
value  '"or  R.  In  the  practical  case  of  a body  containing  many 
dislocations,  R can  in  any  case  only  be  estimated  in  terms  of 
the  average  sjiacing.  of  dislocations. 

Hy  adding  (3),  (5)  and  (12),  we  obtain  for  tlic  total 
volume  of  formation  ['cr  unit  Icng.th  of  dislocation 


where 


6 A 


tot 


c o r e 


a ~ 


+ a...  P.n  (R/h ) + , ] , 

core  M.  ^ outer 

UiP-n(r^yh)  - tn  (r^/h)  , 


”Nh  = ^‘1  ^ 

”outer  " ‘^2  " ^5 


(13) 

(H) 

(15) 

(16) 


In  (13)  we  have  cliosen  h as  the  most  natuiwil  scaling  length 
for  expressing  the  logaritiimic  term.  Tlie  core  related  term 
IS  tlien  g,  iven  hy  the  mode  1 - rcl  axa  t i on  tern  a,  (given  hy  ecys.  1 
and  6)  "corrected"  for  the  effects  of  arbitrary  model  iiarametcrs 
r^  (which  reflects  tiie  relative  positioning  of  atomic  structure 
and  coordinate  system)  and  r^  (tlie  radius  of  the  relaxed  region). 


The  term  reflects  tlic  elastic  properties,  linear  and  non- 

linear, of  the  material,  while  the  last  tenn  includes  the 
effects  of  tlie  outer  boundary  conditions.  A{;ain,  this  term  could 
be  incorporated  in  the  second  by  suitable  choice  of  an  effective 
value  for  R. 

C . Other  flex  Metliods 

The  first  flexible  boundary  metiiod,  wliicli  we  shall  call 
I'lex-S,  was  developed  by  Sinclair^^^.  In  that  method,  atoms 
311  repion  II  and  beyond  are  displaced  from  tlicir  Volterra 
positions,  not  by  Green  functions,  but  by  a tnincated  scries 
of  hiyher  order  elastic  terms  witli  multipolc  s i nyu  1 ^ir  i t i es 
at  the  dislocation  line.  The  i ncomj'a t Lb i 1 i ty  forces  in  region 
II  are  converted  into  generalised  forces  on  tlie  higher  tern 
coefficients,  and  tlicse  forces  with  the  atomic  forces  in  region 
I arc  simultaneously  relaxed.  This  fully  anisotropic  method 
gives  results  in  good  agreement  with.  I'lex  IT.  However,  i'lex-Il 
has  ]iroved  more  efficient  in  computing  tine,  by  a factor  of 
five,  in  tlie  comparisons  of  section  111  and  also  in  trials  by 
Puls  and  iVoo^"*^.  Tlic  oi'iginal  I'lex-S  calculations  used  a matrix 
method  of  relaxation,  which  retiuired  few  force  calculations, 
but  expensive  matrix  nan i jnil at  ions . Later  calculations  using 
tlie  more  p lac  t i call  1 e conjugate  gradients  lelaxation  method 
reqiiiretl  vei'y  many  fence  calculations.  Tlie  latter  relaxation 
method  aiipears  to  respond  jioorly  to  the  particular  type  of 
disequilibrium  which  is  generated  at  intermediate  steps  by  tlic 
ricx-.S  system.  I'lex-S  has,  however,  been  used  by  S inc  1 ‘>i  r ^ ^ 
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in  crack  model  1 ini',  work  and  performed  mucli  more  efficiently 
tlian  in  the  dislocation  case.  The  reason  for  tlie  difference 
is  not  known. 

Another  flexible  boundary  method,  developed  at  about  t!^c 

(2) 

same  time  as  Flex-S,  was  Flex-1  \ Compatibility  between  the 
regions  was  sought  by  calculating  the  tractions  acting  across 
a circular  patli  lying  within  region  I,  representing  these  by 
a Fourier  scries,  and  tJicn  imposing  ailditional  elastic  dis- 
placements on  tlie  outer  region,  corresiioml  ing  to  the  traction 
distribution.  This  method  was  develoju’d  only  in  the  isotrojiic 
approximation.  However,  Teodosiu  et  al^^^^^  liave  recently 
jiroviiied  the  necessary  anisotropic  clastic  solutions,  wliich 
Pctrasch  and  co-workers  liave  used  in  a model  which  is 

essentially  Flex-I  starting  from  non-linear  anisotrojiic  initial 
positions.  However,  we  liave  found  I'lex-I  to  be  much  less 
efficient  tlian  eitlier  lUcx-S  or  Flex-I  1. 

D . Fxtens i on  to  C r a c k s 

For  planar  cracks  witli  straiglit  tijis,  essentially  the  same 
Flcx-11  tcclinitiue  can  be  apjilied.  Tlie  major  difference  is  that 
the  Green  function  ajiproiiriate  for  a cracked  body  must  be  used. 

This  has  been  given  for  general  anisotro;iy  in  the  scxtic  formalism 
by  Sinclair  and  Hirth  and  in  the  isotropic  apjiroximat  ion  by 

('7,0 

Hirth  ct  al.  “ . The  modifications  necessary  because  of  tlie 

non-availability  of  a general  near-neighbour  lattice  Green  function, 
have  already  been  given  at  the  end  of  section  Il.\. 


a 
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III.  R vaults  i 

For  purposes  of  comparison,  we  treat  the  [100]  edge 

dislocation  in  bcc  iion  as  represented  by  the  Johnson-I 
(7) 

potential^  \ This  dislocation  was  already  used  as  the  example  | 

case  for  Flex-S^^^  and  Flcx-I^^^.  ' 

A . I-lex-TI  bfficicncy  ' 

As  an  indication  of  tlic  numerical  efficiency  of  the  metliod, 

we  give  an  example  of  a typical  calculation  on  the  above  ‘ 

dislocation.  Region  I contained  152  atoms,  region  II  a further  i 

78  atoms.  The  desired  accuracy  for  convergence  was  that  | 

- 3 ° - 1 

all  forces  be  smaller  than  10  eVA  . The  conjugate  gradient  ' 

relaxation  method  (CGR)  was  used,  with  a limit  of  50  on  the 

calls  to  the  force  calculating  routine.  The  liistory  of  the 

calculation  was  as  follows.  CGR  (limit  of  50  reached)  on  region 

0-1 

I,  left  forces  in  region  II  of  0.05  eVA  typical;  Green 

functions  applied;  CGR  (50  force  calls):  rcgion-II  forces  0.005 

typical,  0.02  maximum;  Green  function  applied;  CGR(16);  Green 

functions;  CGR (5);  all  forces  in  regions  I and  II  now  less  than 
- 3 ° - 1 

10  ' eVA  ; total  force  calls  = 110.  Tor  com]Kirjson  a run 
of  Flex-S  on  the  same  model  required  426  force  calls  for  tlic 
same  accuracy.  The  original  IUcx-1  calculations  typically 
required  over  1000  force  calculations. 

H . Core  (!on  !'  i gu ra  t i on 

Tlie  dislocation  core  configurations  generated  by  different 
computational  algoritlims  can  be  compared  in  terms  of  the  distances 
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between  atom  pairs  closest  to  the  core  that  are  most  affected 
by  tlic  relaxation.  Table  I presents  a few  selected  interatomic 
distances  for  the  rigid  boundary,  Tlex-S,  Plcx-I  and  Plex-II 
models.  The  atom  pairs  are  those  depicted  in  figure  2.  The 
results  of  the  various  flexible  boundary  treatments  are  in  good 
agreement  with  each  other  and  w'ith  the  large  rigid  boundary 
model.  The  small  differences  between  the  flexible  models  are 
attributable  to  the  varying  termination  criteria  for  the 
relaxation  procedures,  and  to  the  small  residual  errors  at  the 
boundarv,  which  differ  according  to  the  linear  elastic 

f 1 2 251 

approacli  and  region  sliape  in  use.  As  found  in  previous  work,^  ’ ’ ' 

the  effect  of  using  isotropic  rather  than  anisotropic  elasticity 
in  tlie  boundary  conditions  is  not  large;  the  core  configuration  is 
determined  principally  by  the  choice  of  potential. 

Altliougli  not  displayed  in  Table  I,  tlie  dependence  of  core 
atom  sjnacings  on  model  sice,  using  Flex-II,  is  very  small  (<0.05''o) 
down  to  about  100  atoms  in  region  I,  just  as  was  found  witji 

f 1 1 r 2 1 

I'lcx-S''  ‘ and  Flcx-I^‘"  . By  contrast,  witli  rigid  boundary  conditions, 
the  core  atom  spacings  vary  by  several  percent  for  region  I 
in  the  range  90  to  780  atoms 

C.  Dilatation  with  Different  Voltcrra  Field  Rep r es en t a t i ons 

One  of  tlie  aims  of  tlie  present  work  was  to  resolve  a puzzling 
apparent  inconsistency  between  the  Flex-S  and  Flex- I models 
with  resjTcct  to  tlie  dilatational  effect  of  tlie  core.  From  the 
published  values  of  tlie  higlicr-term  coefficients  in  tlie  relaxed 
Flex-S  lioundary,  one  can  calculate  the  area  expansion  of  a 
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r 

circul.-ir  patli  in  the  continuum  rcpjon,  centered  on  the  dislocation. 

f • 

The  result,  for  tlic  cliangc  from  tlie  \'olterra  configuration  to 

[ 

2 

I ll\e  final  relaxed  state,  is  6A  = 0.4Sb  , The  same  quantity  was 

i measured  directly  in  the  flex- I model  from  the  atomic  displacements, 

2 

‘ with  tlie  result  6A  = 0.25b  . In  confirmation,  the  (unpublished) 

coefficients  of  the  higher  terms  in  the  I'lex-I  boundary  fields 

2 

lead  to  tlie  value  6A  = 0.27b  in  tlie  continuum  region. 

In  the  liglit  of  the  analysis  in  section  IIB,  we  now 
recognise  tliat  this  difference  originates  from  the  geometric 
nonlinear  volume  change  associated  witli  the  change  from  perfect 
lattice  to  Volterra  configuration.  Wlicn  Lagrangian  coordinates 
are  used  fas  in  I'lex-I),  the  exact  initial  configuration  actually 
depends  upon  the  angle  at  which  the  plane  of  discontinuity 
in  the  displacement  field  lies.  I'.'itli  bulerian  coordinates,  the 
angle  of  this  branch  cut  is  immaterial,  but  the  configuration 
is  still  sensitive  to  arbitrary  origin  shifts  which  appear, 
from  the  far  field,  to  be  dislocation  climb. 

To  fully  investigate  these  effects,  we  have  used  Flex-II 
to  perform  relaxations  from  four  different  initial  configurations, 
obtained  from  the  anisotropic  Volterra  field  by  the  means 
indicated  in  figure  7> . The  resulting  relaxed  configurations 
were  almost  identical,  as  revealed  by  the  interatomic  spacings 
not  only  near  the  core,  but  also  in  tlie  boundary  region.  However, 
rather  tlian  comparing  results  in  these  terms,  we  display  in  Table 
II  a (juantity  more  directly  related  to  the  question  of  dilatation: 
the  area  of  a local  circuit  enclosing  tlie  dislocation.  The 
circuit  cliosen  was  roughly  square,  passing  tlirough  a sequence 
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of  atom  rows  just  inside  region  1.  The  reference  perfect- 
lattice  area  is  that  for  the  same  number  of  atoms  as  in  the 
circuit,  with  edge  and  corner  atoms  counted  fractionally.  The 
results  of  Table  II  show  firstly  that  tlie  final  configuration  is, 
to  good  approximation,  independent  of  the  path  taken  to  create 
and  relax  the  dislocation.  Secondly,  it  may  be  seen  that 
the  values  quoted  above  for  core  dilatations  in  the  earlier  Flex-S 
and  Flex- I calculations  agree  well  with  the  relaxation  area 
changes  in  the  corresponding  present  results:  compare  case  (d) 
with  the  original  Flex-S  result;  case  (a)  with  Flex-I.  It  is 
important  to  realise  that  the  relative  circuit  area  changes 
cannot  be  associated  witli  the  overall  dilatation  of  the  dislocated 
body,  because 

(i)  relaxation  of  the  outer  surface  of  tl\e  body  leads  to 
an  overall  dilatation  different  from  tliat  of  a circuit 
far  inside  tlie  surface; 

(ii)  geometric  and  clastic  nonl incar i t i es  outside  the  circuit 
contribute  to  tlie  overall  dilatation; 

(iii)  the  circuit  is  not,  and  indeed  could  not  be,  perfectly 
centered  on  the  dislocation,  so  that  even  from  an 
infinitesimal  strain  viewpoint,  the  Volt  err a field  would 
contribute  to  the  area  change. 

We  now  use  the  anal\'sis  of  section  II-B  to  calculate  tltc 
total  formation  of  the  dislocation,  using  the  lUex-11  results 
and  Fuilerian  coordinates.  IVhilc  doing  so,  we  cxiilorc  one  further 
variation  in  tlie  manner  of  creating  initial  positions.  In 


case  (J)  of  Table  II,  atom  positions  in  A were  used  in  the 
1 opa  r i t !nii  functions  of  tlie  Voltcrra  displacement  field.  This 
choice  of  units  is  arbitrary,  but  changing;  the  units  adds  a 
constant  to  the  X2  displacement  function,  which  in  Eulerian 
coordinates  amount  to  shifting  the  initial  perfect  lattice 
relative  to  the  coordinate  origin.  This,  through  altering  Tq 
defined  following  eq.  (3),  alters  tlic  volume  of  the  dislocated 
initial  configuration.  Hence,  to  sliow  that  tlie  final  result 
is  independent  of  this  arbitrary  clioice,  we  made  a further 
calculation  in  which  a constant  additional  X2 ‘displacement  of 
0.1b  was  added  to  the  Voltcrra  field,  equivalent  to  taking  r^  = 

O 

7.293A  in  equation  (3).  After  completing  the  relaxation  in 
each  case,  using  Elex-II,  the  forces  in  region  III  were 
examined  (extending  the  atomic  structvirc  as  necessary]  to  clicck 
the  transition  to  the  distribution  of  forces  predicted  by 
nonlinear  elasticity  tlieory  (eq.  8).  When  moments  of  these 
very  small  forces  were  calculated  as  in  eci.  (1),  taking  the 
sum  out  to  a certain  radius  r,  .the  variation  of  the  sum  with 
Unr  became  linear,  witli  tlie  slope  predicted  from  nonlinear 
elasticity,  just  as  soon  as  the  circle  of  radius  r lay  completely 
outside  region  II,  which  was  square  in  shape.  From  the  intercept 
of  the  graph  against  Unr,  the  effective  value  of  r^  in  eq . (12) 

could  be  found.  Witli  circular  regions  in  use,  r^  could  be 
taken  as  the  inner  radius  of  region  III. 

Tlic  results  which  enter  into  tlie  total  dilatation  are  shown 
in  Table  III.  It  can  be  seen  that  the  effective  shift  of 
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orit^in  by  0.1b  of  tlie  Vollcrra  field  with  I'cspcct  to  the  crystal 
lattice,  leads  to  chani;es  in  the  relaxation  field  (as 
represented  by  moincnts  and  M22  the  dilatation  parameter 

03)  Khicli  almost  exactly  compensate  for  the  difference  of 
O.lb^  in  initial  dilatation,  so  that  the  total  volume  change, 
given  by  and  a^^,,  , is  almost  exactly  the  same.  Tlie  final 
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result,  tlien,  for  t)ie  volume  change  per  unit  length  of  a cylinder 
of  effective  radius  R (sec  Sect.  IIB)  upon  creation  of  a dis- 
location on  the  axis,  is  6A  = [0.37  + 0.070  ln(R/b)]b^.  The 
result  is,  of  course,  dependent  on  the  interatomic  potential 
acting.  I'or  a typical  effective  value  R = lO'^b,  tlie  result  becomes 
6A  = l.Olb^. 

In  order  to  emi)hasize  the  importance  of  the  results  of  this 
section,  we  offer  the  following  physical  description.  The  final 
configuration  of  a given  group  of  atoms  containing  the  dislocation 
is  the  same,  independent  of  the  method  used  to  determine  it. 

For  example,  in  tlie  explored  adjustment  of  the  initial  configuration 
whicii  involved  an  origin  shift  of  0.1b  in  the  climb  direction 
away  from  the  physically  more  obvious  choice  midway  between  two 
atom  planes,  initial  interatomic  forces  were  produced  in  region  I 
such  that  upon  relaxation  tlie  dislocation  "climbed  down"  again. 

This  in  turn  led  to  different  flexible  boundary  relaxations,  so 
that  the  relaxation  volume  change  iSA^^p,  exiircsscd  by  tlie  jiarameter 

7 

03  derived  from  rej;ion-Il  force  di]iole  moments,  was  O.lb"  smaller. 
Ibis  difference  in  tlie  relaxation  volume  was  however  almost 
exactly  compensated  by  that  in  the  non-linear  geometric  volume 
cliange  involveil  in  formation  of  the  initial  configuration.  This 
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componsat  ion  is  expressed  by  ec].  (14),  in  Kliich  tlic  term  in  r^ 
expresses  the  initial  displacement,  in  the  climb  direction,  of 
the  ajijiarent  slip  plane  from  the  coordinate  origin. 

Because  tliis  careful  analysis  is  necessary,  it  may  be  thought 
that  a simpler  approach,  neglecting  only  the  effects  of  far- 
field  clastic  non-linearities,  might  make  use  of  the  area  changes 
of  Burgers  type  circuits.  Indeed,  kc  have  here  shown  that  the 
change,  from  perfect  lattice  to  final  relaxed  state,  of  the 
area  of  such  a circuit  is  independent  of  the  jiatli  taken  in 
creating  and  relaxing  tlic  model.  Allowance  for  outer-surface 
relaxation  would  be  necessary  to  convert  such  circuit  expansions 
into  finite-body  dilatations,  but  a more  serious  difficulty  is 
the  dependence  on  circuit  sice  and  position.  Tlie  above  results 
show  that  the  circuit  would  have  to  be  centered  relative  to  the 
final  dislocation  line  jiosition  and  not  relative  to  the  coordinate 
origin,  which  is  arbitrary.  hven  so  no  choice  for  circviit 
shape  or  position  will  guarantee  correspondence  l)etwecn  tliis 
method  of  estimating  tlie  formation  volume  and  the  more  exact 
analysis.  finally,  it  is  w'orth  emjihasicing  that  in  the  dislocation 
system  investigated  licrc,  the  elastic  nonlinear  term  is 

comparable  to  tlie  core  term,  so  that  the  full  nonlinear  analysis 
is  in  any  case  required, to  give  a realistic  estimate  of  the 
total  ililatation  strength. 


25. 


• Core  I'  ij-l d Rcju~ esc iv^a tie n In’  Li n c J^o jx c P i jioJ  c^s 

A final  jioint  of  interest  is  tlie  description  of  the  core 
generated  field  in  terms  of  a pair  of  orthogonal  line  force 
dipoles  displaced  from  the  dislocation  line.  The  Flex-I 
calculations  found  such  a representation,  witli  the  t^vo  dipoles 
at  the  same  position,  by  a best  fit  to  the  relaxation  displace- 
ments of  atoms  near  the  region  I outer  boundary.  The  present 
calculations  allowed  such  a fit  to  be  olitained  directly  from 
the  dipole  and  quadrupole  moments  of  the  accumulated  Green 
function  source  strengths,  wliich  totally  characterize  the  displace- 
ment field  outside  the  core.  Again,  however,  it  must  be  stressed 
that  tlie  result  refers  only  to  the  relaxat  i on  field,  and  is 
dependent  on  the  choice  of  initial  positions.  Table  IV  compares 
the  Flex- I values  with  the  j^resent  results  using  either  Lagrange 
or  Luler  initial  configurations.  The  rather  poor  agreement 
between  the  Flex-I  and  corresponding,  Flex- 1 I results  jirobably 
reflects  the  different  fitting  criteria:  the  Flex-I  fit  was  near 
to  the  core,  while  the  derivation  from  dipole  and  quadrupole  moments 
rejiresents  a fit  at  "infinity". 


IV.  SUMMARY 

A tcclinique,  Flex-II,  has  been  described  for  ajqilying 
flexible  boundaries  to  an  atomic  model  of  a dislocation,  crack  tip 
or  otlier  line  defect,  so  that  compatibility  is  acliieved  between 
the  core  region  and  tlie  outer  region  described  in  terms  of 
continuum  clastic  solutions.  The  method  has  been  found  to  be 


superior  to  earlier  flexible  boundary  schemes  with  respect  to 
speed  of  computation.  In  the  analysis  of  the  core  contributions 
to  a dislocation  strain  field,  nonlinear  effects  have  been  fully 
allowed  for,  without  tlie  necessity  of  employing  a nonlinear 
elastic  solution  for  tlie  displacements.  This  analysis  has 
enabled  calculation  of  the  overall  dilatation  of  a finite  body 
due  to  a dislocation,  and  earlier  conflicting,  calculations  of 
this  quantity  have  been  corrected  and  reconciled. 
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Tabic  II 

A.'ca  of  Circuit  Around  Pi  si  cat  i o_n  , J'nc  1 os  i n>’  ] 3 4 Atom 
Rows  , \s  i tTr~Tn~l  t~l ; i l~Tos  \ tT o n^s~Ti c i ^cratcwl  ;is  in  Tr]furc~3a~ cT 


Relaxation  by  I'lex-II  with  152  atoms  in  rejjion  I;  Areas,  in  units 
of  b ^ , may  be  cginpared  with  perfect  lattice  value  67.25.  Burgers 
vector  b = 2.8  6.A. 


Case  Area  in  Area  after 

(see  initial  relaxation 

I-  i g . 3 ) c 0 n f i g u r a t i 0 n 


Change  during 
relaxation 


67.4257 
67  . ; 
67.4262 
67 . 2537 


67, 

, 678  7 

67  , 

,6776 

67. 

, 6818 

67. 

,6695 

0.2529 

0.4054 

0.2555 

0.4158 


Table  III 

Contributions  to  Total  Volume  of  Tormation  of  Dislocation 


luo  cases  (sec  text)  using  liulerian  Volterra  solution  for  initial 
positions. 


Case 

Fig.  3(d) 

Fig.  3(d)  i t h 
0.1b  added  to  U2 

"1 

-0 . 0503 

-0.0503 

0 

aiinCi-Q/A) 

0 

-0.1 

Mj^/evA-^ 

6.876 

5.884 

M22/cVX‘^ 

2.725 

2.  238 

0.6521 

0. 5517 

"4 

0.1204 

0.1204 

>o 

19.1 

19.2 

“core 

o 

o 

0.3696 

“NL 

0.0701 

0.0701 

1 
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Table  IV 

Dcscr  ijit  ion  of  Relaxation  Pi  spl  accnent  Field  in 
^ ^ ’ 0 f e Dipoles 

O •»  O 

nipolc  streiii^tl'S  are  in  eVA'*^,  positions  in  A.  All  dipoles  have 
coordinate  xj  = 0. 


Case 

Mil 

X2  position 

M22  ^2  position 

of  M2 2 

F 1 e X - 1 as 
fig.  3(a) 

9.4 

-5.2 

2.8  -5.2 

Flex- 1 1 as 
fig.  3(a) 

5.0 

-6.1 

0 

OC 

1 

FI  ex  - 1 1 as 
fig.  -3(d) 

6.9 

-2.5 

2.7  -1.0 

Fif,ure 

ri^:.  1 

Fig.  2 
Fig.  3 


■ a p 1 1 0 n s 

Concolric  arrangement  of  regions  for  Flcx-II  model. 

The  Shape  need  not.  he  circular. 

Relaxed  configuration  of  the  core  of  a [100]  edge 
dislocation  in  hcc  iron.  Two  (001)  layers  are  distinguished. 
Four  different  methods  of  obtaining  initial  dislocated 
configuration  from  the  Volterra  solution.  (a),  (b) 
and  (c)  use  I.agrangian  coordinates,  with  the  plane  of 
discontinuity  (a)  on  tlie  positive  axis,  witli  a lialf 
plane  of  atoms  inserted;  (b)  on  the  x^  axis,  so  that  slip 
is  the  formation  im'occss;  (c)  on  tlie  negative  X2  axis, 
with  a lialf  jilane  of  atoms  removed.  In  (d)  , liulcrian 
coordinates  arc  used,  and  the  location  of  the  plane 
of  discontinuity  is  immaterial. 


